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ABSTRACT 


The  calculation  of  non-lifting  potential  flow  about  arbitrary  three 
dimensional  bodies  is  examined  in  detail  with  specific  interest  in  the  XYZ 
Potential  Flow  program  developed  by  the  David  W.  Taylor  Naval  Ship 
Research  and  Development  Center.  The  program  uses  a  surface  singularity 
distribution  to  solve  the  Neumann  boundary  value  problem  by  means  of  a 
source  panel  method  assuming  a  flat  element  with  a  constant  source 
density  over  the  area  of  the  element.  Boundary  conditions  are  applied  at 
control  points  on  the  elements  producing  a  system  of  linear  equations  for 
the  source  density.  When  the  source  density  is  known,  velocities  and 
pressure  coefficients  may  be  calculated. 

The  main  purpose  of  this  paper  is  to  present  the  details  of  the 
approximation  of  an  arbitrary  three  dimensional  body  using  quadrilateral 
elements,  and  to  provide  a  detailed  derivation  of  the  exact  source  panel 
integrations  in  order  to  gain  insight  for  future  research  at  Texas  A&M 
University.  A  variation  of  the  Hess  method  of  surface  discretization 
using  quadrilateral  source  panels  is  described  in  detail  as  it  is  used  in  the 
XYZ  Potential  Flow  program.  The  exact  source  panel  integrations  are 
derived  in  detail. 

A  general  discussion  of  other  aspects  of  the  program  is  included. 
Velocities  and  pressure  coefficients  for  flow  about  a  triaxial  ellipsoid 
are  calculated  using  the  XYZ  Potential  Flow  Program  solution,  and  the 
results  are  compared  with  the  analytical  solution  and  the  Hess  Program 


solution. 
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This  paper  examines  two  aspects  of  the  development  of  the  XYZ 
Potential  Flow  Program  (hereafter  referred  to  as  the  XYZPF  Program),  a 
FORTRAN  program  which  uses  a  source  panel  method  to  approximate 
solutions  to  steady  potential  flow  problems  about  arbitrary  three 
dimensional  bodies.  The  aspects  examined  in  detail  are  (t)  the 
description  of  the  details  of  the  approximation  of  an  arbitrary 
three-dimensional  body  using  quadrilateral  elements,  and  (2)  a  detailed 
derivation  of  the  exact  source  panel  integrations. 

The  XYZPF  Program  was  developed  specifically  for  applications  in 
numerical  ship  modelling  and  hydrodynamics  studies  at  the  David  W. 
Taylor  Naval  Ship  Research  and  Development  Center  (NSRDC)  in  Bethesda, 
Maryland.  The  format  of  the  program  is  based  on  the  work  of  Hess  and 
Smith  (1962)  in  the  numerical  calculation  of  non-lifting  potential  flow. 
A  similar  program  is  maintained  by  the  Aerodynamics  Division  of  the 
McDonnell -Douglas  Corporation,  referred  to  in  this  paper  as  the  "Hess 
program."  The  XYZPF  Program  is  a  modification  of  what  has  come  to  be 
known  generally  as  the  Hess  Method.  The  most  significant  modifications 
are  improvements  to  the  method  of  solving  the  matrix  equation  for  the 
source  density,  and  greater  flexibility  in  the  input  options  (Dawson  and 
Dean  1972). 

Though  potential  flow  is  a  product  of  mathematics,  and  is  never 
found  in  a  real  fluid,  the  results  of  potential  flow  calculations  provide 
usable  information  for  flow  regions  external  to  a  thin  boundary  layer, 


with  little  or  no  boundary  layer  separation.  For  such  flow  fields,  the 
region  outside  the  boundary  layer  may  be  considered  to  be  effectively 
inviscid,  and  may  be  closely  approximated  by  potential  flow  models. 

Small  viscous  effects  can  be  accounted  for  by  "thickening"  the  body  by  the 
appropriate  displacement  thickness.  Displacement  thickness  accounts  for 
the  region  of  retarded  fluid  flow  in  the  boundary  layer  inversely 
proportional  to  the  square  of  the  free  stream  velocity.  Downstream  of 
the  point  of  boundary  layer  separation,  the  potential  flow  model  no 
longer  applies. 

Prior  to  the  development  of  numerical  methods,  analytical 
solutions  were  generally  restricted  to  simple  analytical  shapes  (Kellogg 
1929).  The  need  to  solve  boundary  value  problems  for  arbitrary 
boundaries  in  continuum  mechanics  has  fostered  the  development  of 
numerical  approximations  to  the  integral  equation  expressions.  While 
the  integration  methods  have  been  well  known  for  quite  sometime,  only 
since  the  advent  of  high  speed  computers  have  many  of  the  problems 
been  practical  to  solve  by  numerical  methods.  Among  the  numerical 
methods  being  used  are  finite  differences,  finite  elements,  and  the 
boundary  element  method. 

"Finite  differences”  and  "finite  elements"  are  numerical  methods 
which  satisfy  the  boundary  conditions,  and  then  approximate  the 
solution  to  the  governing  equation  in  the  fluid  domain.  These  methods 
discretize  the  domain  into  a  network  of  elements  or  cells. 
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Another  approach  is  what  is  now  known  as  the  "Boundary  Element 
Method,"  in  which  the  governing  equation  is  exactly  satisfied  in  the 
domain,  and  the  boundary  conditions  are  applied  through  a  boundary 
discretization  method.  The  boundary  value  problem  is  reformulated  as  a 
boundary  integral  equation  which  is  then  discretized  by  subdividing  the 
boundary  into  a  finite  number  of  surface  elements.  Each  element  is 
represented  by  an  analytical  function,  and  the  source  density  function  is 
integrated  over  the  surface  of  each  element.  Two  factors  governing  the 
accuracy  of  the  boundary  element  method  are  the  boundary  discretization 
method  and  the  source  panel  integration.  These  two  factors  are  examined 
in  detail  in  this  report,  as  a  detailed  derivation  of  the  exact  source  panel 
integration,  including  the  development  of  the  source  panel  geometry,  has 
not  previously  appeared  in  literature. 

The  difference  between  the  domain  methods  and  the  boundary 
methods  is  significant.  The  domain  methods  discretize  the  domain,  while 
the  boundary  methods  discretize  the  boundary.  Thus,  the  boundary  method 
reduces  the  dimension  of  the  problem  by  one,  as  depicted  in  figure  (1).  In 
the  application  of  the  XYZPF  Program,  the  problem  is  reduced  from  a 
three-dimensional  problem  in  the  domain  to  a  two-dimensional  problem 
on  the  boundaries.  This  method  is  well  suited  to  problems  in  which  the 
limits  of  the  domain  are  infinite  or  difficult  to  define,  in  that  the 
problem  is  applied  to  the  boundary  rather  than  the  domain. 
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FINITE  ELEMENT  DISCRETIZATION 


BOUNDARY  ELEMENT  DISCRETIZATION 


Ftf«r*  1 .  Discretization  Methods 

Just  as  there  are  many  variations  of  domain  methods,  there  are 
also  a  variety  of  boundary  methods.  In  general,  they  can  be  classified  as 
"indirect"  or  "direct"  formulations.  The  "indirect"  method  assumes  a 
continuous  source  distribution  over  the  surface  of  the  body,  and  a 
solution  which  satisfies  both  the  governing  equation  in  the  domain,  and 
the  boundary  conditions  on  the  body  surface.  The  result  is  an  integral 
equation  on  the  boundary  which  has  the  surface  source  density  function  as 
its  unknown.  By  enforcing  the  boundary  conditions  at  control  points  on 
the  surface,  a  system  of  equations  is  produced  by  which  the  source 
density  may  be  determined. 

The  "direct"  method  solves  the  velocity  potential  function  through 
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an  application  of  Green's  Second  Identity  requiring  the  solution  of  a 
source  distribution  and  a  dipole  distribution  on  the  boundary.  The  direct 
method  has  more  physical  significance  to  the  general  boundary  value 
problem,  and  more  versatility  in  its  application  as  it  can  be  applied  to 
Neumam  problems,  Dirichlet  problems,  or  mixed  boundary  value  problems 
(Brebbia  1984). 

The  simplicity  and  accuracy  of  the  indirect  method  has  made  it 
attractive  for  many  applications  The  source  panel  method  is  an 
application  of  the  indirect  formulation  of  the  boundary  element  method 
to  the  Neumann  type  of  potential  flow  problem,  for  which  the  normal 
derivative  of  the  potential  function  is  prescribed  on  the  boundary. 

1.1  OBJECTIVES 

The  purpose  of  this  paper  is  (1)  to  describe  the  details  of  the 
approximation  of  an  arbitrary  three-dimensional  body  using  quadrilateral 
elements,  and  (2)  to  provide  a  detailed  derivation  of  the  exact  source 
panel  integrations  for  use  in  future  investigations  at  Texas  A&fi 
University  using  panels  of  higher  order  geometries  and  source  density 
functions.  This  paper  is  not  intended  to  be  a  user’s  manual,  though  a 
general  discussion  of  other  aspects  of  the  program  is  also  included. 
NSRDC  Report  3892  (Dawson  and  Dean  1972)  is  a  summary  of  the  XYZPF 
Program  for  those  strictly  interested  in  its  use. 


The  foundations  of  the  boundary  element  method  were  laid  early  in 
this  century  beginning  with  Fredholm  in  1903  when  he  established  the 
existence  of  solutions  to  the  Neumann  problem  through  a  reconstruction 
of  the  problem  using  a  discretized  boundary  (Kellogg  1929).  The  solution 
was  determined  to  be  the  potential  of  a  simple  source  distribution  on  a 
boundary  with  a  continuous  normal  derivative  for  an  infinite  domain. 

Later  works  by  Kellogg  (1929)  in  potential  theory  demonstrated  the 
application  of  the  boundary  integral  equation  method  in  electrostatics, 
heat  transfer,  flow  in  porous  media,  and  fluid  flow  problems,  but 
development  was  limited  by  the  difficulty  of  obtaining  analytical 
solutions.  No  significant  advances  were  made  until  interest  inboundary 
integral  equation  methods  was  revitalized  with  the  advent  of  high  speed 
electronic  computers.  Investigators  were  then  able  to  discretize  the 
boundaries  and  solve  the  integral  equations  numerically.  This  method  of 
solution  became  known  as  the  boundary  element  method.  Early 
development  of  such  numerical  methods  was  pioneered  by  Hess  and  Smith 
(1962)  and  Jaswon  and  Symm  (1963).  Hess  and  Smith  dealt  primarily  with 
the  indirect  formulation  eventually  leading  to  a  solution  for  the  three 
dimensional  problem  as  described  in  this  paper.  In  a  parallel  work, 
Jaswon  and  Symm  developed  a  direct  formulation  approach  to  the  two 
dimensional  problem.  TheXYZPF  Program  is  based  primarily  on  the  work 
of  Hess  and  Smith.  Hess  has  since  developed  a  higher  order  panel  method 
(Hess  1979)  and  Lefebvre  modified  the  XYZPF  Program  for  calculating 
velocity  potentials  for  five  degrees  of  freedom  (Lefebvre  1982). 
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3.1  THE  POTENTIAL  FLOW  PROBLEM  IN  THREE  DIMENSIONS 

The  governing  equation  for  idea!  (incompressible,  inviscid, 
irrotationa!)  flow  is  Laplace's  equation: 

V2*  =  0  0) 

where  §>  is  the  velocity  potential,  and  V2  is  the  Laplacian  operator.  The 
XYZPF  Program  deals  with  steady,  uniform  flow  of  an  ideal  fluid  about  an 
arbitrary  three  dimensional  body.  The  velocity  components  at  any  point 
within  the  flow  field  may  be  obtained  from  the  negative  gradient  of  the 
velocity  potential,  i.  e. 


V  =  -V* 


ax  ~  ay  J  az 


The  freestream  flow  V~  is  defined  as  a  uniform  stream  of  unit 


magnitude. 


o  n  O' 

V  +  v4 ,,  +  V"  =1 
oox  °°y  wz 


The  key  to  the  boundary  element  method  is  the  Divergence  Theorem 
(Green's  Theorem)  which  relates  a  volume  integral  to  an  equivalent 
surface  integral  reducing  the  three-dimensional  problem  to  a 


two-dimensional  one.  The  expression  for  Green's  second  identity  is 
(Lamb  1924): 


-  wV2$)  dft  =  JJ(w|^  -  $  )  dr  (4) 

in  which  Q  represents  the  integration  over  the  three  dimensional  domain, 
and  T  represents  integration  over  the  two  dimensional  boundary.  The 
partial  derivatives  are  taken  with  respect  to  the  outward  normal,  n.  The 
weighting  function,  w,  is  usually  chosen  to  be  the  fundamental  solution 
for  three  dimensions,  w  =  1/(4rrr),  where  r  is  the  distance  from  the 
source  to  an  arbitrary  point  on  the  boundary. 


Fifvre  2.  PotentUl  Flov  in  Three  Dimensions 


Consider  an  arbitrary  three-dimensional  body  with  surface  T, 
having  an  equation  of  the  form  F(x,  y,  z)  =  0  where  x,  y,  z  are  Cartesian 
coordinates  of  the  global  reference  system  as  shown  in  Figure  (2).  The 
unit  outward  normal,  n,  at  any  point  on  the  surface  is  given  by  the 
gradient  of  the  function  describing  the  surface  divided  by  the  magnitude 


of  the  gradient,  i.e. 


n 


±VF 

VF 


(5) 


where  the  sign  of  the  unit  normal  vector  is  chosen  to  ensure  that  the 
vector  is  an  outward  normal.  The  potential  function  $  describing  the 
flow  field  must  meet  the  following  boundary  conditions: 

a.  V2*  =  0  (Laplace's  Equation)  (6) 

b.  For  an  impermeable  boundary,  the  velocity  normal  to  the  surface 
must  be  zero  relative  to  the  boundary  (the  Neumann  boundary 
condition): 


c.  The  velocity  potential  approachesthe  freestream  velocity 
potential  as  the  distance  from  the  body  goes  to  infinity: 

*  *«,  as  |r  |  -►  co  (8) 

The  total  potential  at  any  point  In  the  domain  Is  composed  of  the 
freestream  potential  and  the  disturbance  potential  due  to  the  body, 

♦  =  t„  *  t  (9) 
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The  disturbance  potential,  <p,  satisfies  the  following  boundary 
conditions: 

a.  V2<p  =  0  (10) 

b.  From  equation  (7),  the  velocity  normal  to  the  boundary  due  to 
the  disturbance  and  due  to  the  onset  flow  must  be  of  equal 
magnitude,  but  opposite  sign.  Then  from  equation  (9) 

=  n(p)  •  Y*,  (ID 

V  On  Jr 

Note  that  the  normal  vector  is  a  function  of  position  on  the  surface 
of  the  body 

c.  The  disturbance  potential  approaches  zero  as  the  distance  from 
the  body  goes  to  infinity,  i.  e. 

9  — ►  0  as  |  r  |  — ►  oo  (12) 

3.2  MATHEMATICAL  MODEL 

The  disturbance  potential  of  the  body  may  be  represented  by  a 
distribution  of  a  source  density  function  o  over  the  body  surface.  The 
potential  at  an  arbitrary  point  P(x,  y,  z)  due  to  the  surface  potential  is 
(Kellogg  1929): 


.  *  ■»  ^  v* ;■*  .  ■  .-*¥  ■  ^  TO 


f  (X,  y,  2)  = 


(13) 


wher.  q  is  the  integration  point  on  the  surface,  and 


r(P,q)  = 


(x  -  x0r  ♦  (y  -  y0)  *  (z  -  z0)' 


is  the  distance  from  the  field  point  P  to  the  integration  point  q. 


The  source  density  distribution  function  must  satisfy  the  boundary 
conditions  for  the  disturbance  potential.  Boundary  conditions  (10)  and 
(12)  are  automatically  satisfied  by  the  form  of  the  integrand.  However, 
equation  (11),  the  velocity  normal  to  the  boundary,  combined  with  the 
Neumann  boundary  condition,  equation  (7),  is  the  key  to  solving  the 
boundary  integral  problem. 


The  integrand  becomes  singular  as  the  surface  of  the  body  is 
approached.!,  e.  |r|  goes  to  zero.  The  singularity  represents  the  local 
fluid  flux  normal  to  the  boundary  due  to  the  local  source  density.  The 
principal  value  of  the  singularity  is  -2Tro(p),  determined  through  a 
limiting  processof  the  Gauss  Flux  Theorem  (Kellogg  1929).  The  point  p 
represents  a  field  point  which  lies  on  the  boundary.  The  integral 
expression  is  now  composed of  the  contribution  of  the  local  source 
density  and  the  contribution  of  the  source  density  function  over  the 
remainder  of  the  body  surface.  Solving  for  the  velocity  norma!  to  the 
surface  yields  the  following  expression; 


n 


(14) 


¥)r-  -»«'»•  /;» 


dr 


From  equation  (11),  this  expression  becomes; 


2tto(p) 


-n 


o(q) 


r(p,q) 


dr  =  -  n  (p)  ■  V. 


00 


(15) 


This  equation  is  a  two  dimensional  Fredholm  integral  equation  of  the 
second  kind,  which  ensures  a  unique  solution,  and  that  the  diagonal 
elements  of  the  system  matrix  will  be  dominant,  each  having  a  value  of 
2tt  (Kellogg  1929).  Once  equation  (15)  has  been  solved  f or  the  source 
density  o,  the  velocity  components  at  any  point  of  the  flow  field  may  be 
obtained  by  differentiating  the  disturbance  potential  function  (13)  with 
respect  to  the  coordinate  directions  and  adding  the  components  of  the 

freestream  flow,  VM. 


V(x,  y,  z)  =  Voo  - 


5U  (16) 


The  body  shape  does  not  have  to  be  slender,  axisymmetric,  or 
simply  connected,  allowing  for  analysis  of  interior  flow  and  a  wide  range 
of  applications  of  the  method.  The  only  restriction  imposed  on  the  form 
of  the  body  is  that  it  must  have  a  continuous  normal  vector. 
Discontinuities  in  the  right  hand  side  of  equation  (15)  will  produce 
unwanted  singularities.  Thus,  in  the  process  of  approximating  a  body 
which  has  distinct  corners,  where  there  is  clearly  a  discontinuity  in  the 
normal  vector,  the  corner  must  be  replaced  by  a  surface  with  some  finite 


curvature.  However,  application  of  this  method  has  shown  that  the  flow 
calculations  give  correct  results  for  convex  corners,  while  unrounded 
concavecornersmay  or  may  not  producesignificant  error,  depending  on 
the  angle  produced  by  the  corner  (Hess  and  Smith  1962). 


Because  of  the  method  of  approximation,  the  calculation  of  flow 
velocities  on  the  body  surface  are  restricted  to  the  points  at  which  the 
boundary  conditions  were  applied.  Velocities  at  points  other  than  those 
must  be  obtained  by  interpolation.  Direct  calculation  of  velocities  at  the 
edge  of  an  element  yields  infinite  velocities. 

With  the  solution  of  the  system  of  linear  equations  for  the  source 
densities,  the  flow  velocities  at  any  point  in  the  domain  may  be  obtained 
from  equation  (16),  and  pressure  coefficients  are  then  computed  from  the 
velocities  using  a  form  of  the  Bernoulli  equation: 


p«)  =  y  *  tM2  *  07) 

where  P(t)  is  a  constant  independent  of  position.  In  the  XYZPF  Program, 
the  flow  is  steady.  Therefore,  equation  (17)  can  be  reduced  to 


P 


constant 


(18) 


and  the  pressure  field  can  be  expressed  in  terms  of  a  dimensionless 
pressure  coefficient  Cp  as: 


P  -  p, 


oo 


tPV, 


00 


(19) 
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where  p*,  is  the  static  pressure  at  infinity. 


3.3  NUMERICAL  MODEL 


In  order  to  represent  the  surface  of  a  body  in  the  domain 
mathematically,  the  body  may  be  described  by  analytical  expressions 
which  may  provide  an  exact  representation  of  the  surface.  However,  the 
types  of  bodies  which  can  be  adequately  described  by  such  methods  are 
severely  limited.  Another  way  to  represent  the  body  is  to  use  a  large 
number  of  analytical  expressions,  each  describing  only  a  small  portion  of 
the  body.  Hess  and  Smith  (1962)  suggested  the  use  of  an  assembly  of  flat 
quadrilateral  elements  to  model  the  actual  surface  of  the  body,  as  shown 
in  Figure  (3).  Each  quadrilateral  approximates  a  region  of  the  surface 
described  by  points  which  lie  on  the  actual  surface  of  the  body.  As  planar 
elements,  these  quadrilaterals  are  clearly  analytical,  and  when  carefully 
constructed,  the  elements  can  approximate  arbitrary  three  dimensional 
body  surfaces  without  restriction. 


TheXYZPF  Program  uses  the  discretization  procedure  described  by 
Hess  and  Smith  (1962)  with  some  minor  modifications.  The  three 
dimensional  body  surface  may  be  described  using  a  large  number  of  plane 
quadrilateral  elements,  each  assumed  to  have  a  constant  source  density 
over  the  area  of  each  element.  Regions  of  the  body  requiring  higher 
resolution  for  sharp  curvature  or  anticipated  velocity  gradients  will 
require  a  higher  concentration  of  elements. 

Because  the  plane  quadrilateral  elements  cannot  fit  edge  to  edge  on 
a  rounded  surface,  small  gaps  in  the  panel  approximation  contribute  to  the 
error  of  the  approximation.  However,  the  error  due  to  the  gaps  is 
negligible  when  compared  with  the  error  of  the  basic  model,  that  is,  using 
flat  panels  to  approximate  a  curved  surface.  Triangular  elements  have 
been  suggested  in  an  attempt  to  eliminate  the  gaps  (Levy  1959),  but  the 
increased  accuracy  is  so  small  that  it  may  not  justify  the  additional  work 
of  organizing  the  triangular  elements  in  lieu  of  the  simpler  quadrilaterals 
(Hess  and  Smith  1966).  The  method  presented  is  valid  for  an  polygonal 
element  with  any  number  of  sides. 

Equation  (15)  can  now  be  decomposed  into  a  summation  of 
integrals,  each  representing  the  contribution  of  one  element  of  the  body 
surface.  The  unknown  source  density  can  be  taken  outside  the  integral, 
since  it  is  assumed  to  be  a  constant  over  each  element.  The  integration  is 
performed  over  the  area  of  the  source  element,  and  the  boundary 
condition  equation  (11)  is  then  enforced  at  a  single  point  p  in  each 
remaining  element.  By  performing  this  operational  each  element  of  the 


surface,  a  system  of  linear  equations  is  generated  which  is  equal  in 
number  to  the  number  of  surface  elements  and  the  number  of  unknown 
source  densities.  Equation  (15)  can  now  be  approximated  by  the  matrix 
equation  (Dawson  and  Dean  1972): 

a  =  y  a  c  •  +  v  (20) 

.i 

where 
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It  is  important  to  note  that  the  influence  coefficients  Cjj  and  Cjj  are 

functions  of  geometry  only,  and  once  computed,  need  not  be  recomputed 
for  analysis  of  several  different  flows.  From  the  solution  of  equation 
(20)  on  the  discretized  surface,  equation  (13)  may  be  applied  at  any  point 
in  the  domain.  Then,  the  velocity  at  an  arbitrary  field  point  P(x,  y,  z)  in 
the  domain  may  be  determined  from  equation  (16).  With  the  velocity 
known,  the  pressure  coefficient  is  determined  from  equation  (19). 


The  XYZPF  Program  is  actually  composed  of  seven  independent 
programs,  referred  to  as  sections  PF!  through  PF7,  each  of  which  builds 
on  data  generated  from  a  previous  section.  This  type  of  organization 
allows  the  user  the  flexibility  of  rerunning  portions  of  the  program  using 
different  flow  parameters  without  having  to  go  through  the  time 
consuming  process  of  recalculating  the  influence  coefficient  matrix, 
which  is  dependent  only  on  the  geometry  of  the  body.  While  the  NSRDC 
program  is  very  similar  to  the  Hess  program,  there  are  also  some 
significant  differences.  The  following  list  of  differences  is  taken  from 
N5RDC  Report  3892  (Dawson  and  Dean  1972): 

(1)  The  input  to  XYZ  PF  is  arranged  to  facilitate  the  preparation  of 
input  for  a  series  of  problems  in  which  only  one  part  of  the  body  is 
changed.  Also,  a  number  of  checks  are  made  on  the  input  to  help  detect 
errors. 

(2)  An  option  was  added  for  the  recomputation  of  the  source 
density  and  velocities  for  only  part  of  the  body  when  only  small  changes 
are  made.  This  option  also  provides  for  the  use  of  the  solution  of  one 
problem  as  an  initial  guess  for  the  solution  of  another  problem. 

(3)  The  matrix  of  influence  coefficients  is  computed  column  by 
column  instead  of  row  by  row.  This  column  arrangement  was  used  for  the 
original  LARC  computer  version  because  it  required  much  less  high  speed 
memory.  The  computation  is  also  about  10 &  faster  this  way  than  with 
the  row-by-row  arrangement. 

(4)  a  simultaneous  displacement  iterative  scheme  is  used  to  solve 
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the  matrix  coefficient  for  the  source  density.  The  scheme  is  slower  than 
the  successive  displacement  (Gauss-Seidel)  scheme  used  in  the  Hess 
program,  but  it  can  be  carried  out  using  the  matrix  column  by  column 
instead  of  row  by  row. 

(5)  When  possible,  an  extrapolation  procedure  is  used  to  reduce  the 
number  of  iterations  required  for  convergence.  One  such  method  is  the 
Richardson  extrapolation. 

The  methods  used  in  the  XYZPF  Program  will  be  discussed  in  detail 
in  the  following  sections. 


5.0  DETAILS  OF  THE  SURFACE  APPROXIMATION 

5.1  PREPARATION  OF  THE  INPUT 

Section  PFI  is  set  up  to  read  and  store  the  input  data,  and  to 
examine  the  cornerpointsof  the  quadrilaterals  to  detect  obvious  errors  in 
the  input.  Because  of  the  number  of  points  which  may  have  to  be  entered 
for  a  complex  geometry,  the  user  input  is  a  major  source  of  program 
error,  and  this  first  look  for  input  errors  will  save  lot  of  run  time  in  the 
program  as  a  whole.  If  Section  PFI  detects  major  errors  in  the  input,  the 
program  will  not  continue  with  the  calculation  of  the  coefficient  matrix, 
but  will  stop  and  identify  the  grid  location  of  the  error.  Minor  errors  may 
not  cause  the  program  to  stop,  but  will  be  noted  in  the  output. 

One  of  the  major  advantages  to  this  program  is  in  the  organization 
of  the  input  data.  The  surface  is  input  in  sections  so  that  small  portions 
of  the  input  geometry  may  be  changed  without  having  to  recalculated  the 
points  for  the  entire  surface.  The  program  also  takes  advantage  of 
symmetry  to  minimize  the  input  effort.  Only  the  portion  of  the  body 
which  has  no  redundancy  needs  to  be  entered  point  by  point.  The 
remainder  of  the  body  is  reflected  across  the  planes  of  symmetry  by  the 
program  to  complete  the  surface  representation. 

The  surface  is  represented  by  a  set  of  points  in  three-dimensional 
space  which  lie  on  the  actual  surface,  and  which  will  later  be  used  to 
define  the  plane  quadrilateral  source  elements.  These  points  are  defined 
in  the  global  reference  system.  The  points  on  the  surface  should  be 
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selected  in  such  a  way  as  to  provide  an  accurate  representation  of  the 
surface  with  the  fewest  number  of  points  possible.  Portions  of  the 
surface  which  are  highly  curved  will  require  a  larger  number  of  points  to 
provide  adequate  resolution.  Additionally,  portions  of  the  surface  in 
which  the  flow  field  is  expected  to  change  rapidly  will  require  a  large 
number  of  points  to  accurately  determine  the  flow  field  in  that  region. 
Some  familiarity  with  fluid  dynamics  will  provide  a  somewhat  intuitive 
approach  to  properly  distributing  the  elements.  Elements  should  change 
gradually  in  size  from  areas  of  high  concentration  to  those  of  just  a  few 
large  elements,  changing  no  more  than  50  percent  in  size  between 
adjacent  elements  (Hess  and  Smith  1966).  The  accuracy  is  only  as  good  as 
that  provided  by  the  largest  element  in  a  particular  area.  The  use  of 
quadrilateral  elements  facilitates  the  use  of  known  analytical  equations 
and  body  contours  to  determine  the  input  points. 

/ - - - , 


Ftf«re  4.  The  3D  quadrfUtaral  element  In  global  coordinates 

For  the  purposes  of  this  program,  the  body  surface  is  represented 
by  a  large  number  of  plane  quadrilateral  elements  as  shown  in  figure  (3), 
each  of  which  is  assumed  to  have  a  constant  value  of  source  density  over 
the  area  of  the  element.  Each  element  is  defined  by  four  input  points 
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which  lie  on  the  actual  surface  as  shown  in  figure  (4).  Since  each  input 
point  can  be  used  as  a  corner  for  up  to  four  elements,  there  is  no  need  to 
enter  the  same  point  four  separate  times.  The  input  points  are  organized 
in  groups  of  four  to  form  the  quadrilateral  element,  and  each  point  may 
also  be  associated  with  adjacent  quadrilaterals.  This  is  accomplished 
through  the  use  of  a  two  dimensional  coordinate  system  in  which  the  user 
assigns  a  pair  of  integers,  m  and  n,  to  each  point  which  identifies  the 
"row"  and  the  "column”  in  which  it  lies.  A  column  of  points  will  be  given 
a  commonvalue  of  n,  and  each  point  in  that  column  will  have  a  unique 
value  of  m  corresponding  to  the  row  in  which  the  point  lies.  The 
orientation  of  these  "coordinate"  integers  determines  the  direction  of  the 
outward  normal  for  each  element.  Looking  from  the  flow  field  toward 
the  section  of  elements,  if  the  values  of  m  are  increasing  upward,  the 
values  of  n  must  increase  to  the  right.  Increasing  m  and  n  can  point  in  any 
direction  with  respect  to  the  global  reference  system.  In  fact,  the 
orientation  can  change  from  one  section  to  another.  However,  any  other 
relationship  between  m  and  n  will  produce  an  incorrect  normal  vector. 
Once  assigned,  the  values  of  m  and  n  also  serve  to  identify  the  element 
for  which  the  corresponding  point  is  a  corner.  The  four  points  which 
form  a  quadrilateral  element  are  two  points  of  one  column,  or  n  line,  with 
consecutive  m  numbers,  and  two  points  of  the  next  higher  n  line  with  the 
same  m  numbers  as  the  previous  two  points.  Thus,  the  element  m,  n  is 
composedof  the  points  identified  by  (m,  n),  (m+1,  n),  (m,  n+1),  and  (rrH, 
n+1). 
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Ftf«rt  5.  Quadrilateral  index  numbers 

Each  section  of  the  body  surface  is  formed  by  specifying  a  set  of 
cornerpoints  corresponding  to  the  m,  n  pairs  for  all  of  the  quadrilaterals 
of  the  section.  The  user  will  sequentially  assign  an  m  number  to  the 
points  for  each  n  line,  and  also  number  the  n  lines  for  the  section  points 
entered.  The  first  point  in  each  n  line  will  always  have  m  =  1.  The  n  lines 
are  also  numbered  sequentially,  but  the  value  of  n  is  not  reset  for  each 
new  section.  The  sequence  of  n  numbers  runs  through  all  the  sections  as 
shown  in  figure  (5).  Points  on  a  particular  row  or  column  do  not  have  to 
be  strictly  colinear.  By  forming  nearly  triangular  elements,  a  rounded 
planform  can  be  approximated  without  conflicting  with  the  numbering 
convention,  as  shown  in  figure  (6). 
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Ftgor*  6.  Approximating  a  thin  region  vtth  rounded  planform  (Hess  &  Smith  1 962) 


By  entering  data  in  sections,  small  changes  in  geometry  can  be 
performed  without  having  to  reenter  all  the  points  associated  with  the 
body  surface.  This  feature  is  unique  to  the  XYZPF  Program,  and  offers  a 
great  deal  of  flexibility  in  design  work.  However,  with  the  added 
flexibility  comes  more  restrictions  on  both  the  input  of  the  original 
geometry  and  on  any  modifications.  There  are  four  important 
restrictions  on  the  input  which  are  required  to  provide  quadrilateral 
elements  in  groups  of  four  to  facilitate  geometry  calculations  (Dawson 
and  Dean  1972): 


(1)  There  must  be  an  even  number  of  elements  in  both  the  m  and  n 
directions  in  each  section  of  the  body. 


1(2)  The  common  corner  point  of  a  group  of  four  elements  must  not 
coincide  with  any  other  corner  point.  The  sides  between  the  elements 
!.*  serve  to  define  the  local  coordinate  system,  and  serve  as  the  axis  of 

rotation  when  the  surface  is  flattened  for  numerical  differentiation  of 
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the  velocity  potential. 


(3)  Each  set  of  four  elements  must  have  at  least  seven  distinct 
corner  points  to  allow  the  elements  to  more  closely  conform  to  a  curved 
surface.  This  also  allows  for  convergenceof  N-tines  orh-lines  as  may 
occur,  for  example,  at  the  leading  edge  of  an  ellipsoid.  Thus,  only  two  of 
the  four  quadrilateral  may  degenerate  into  triangles  by  having  two  of 
their  corner  points  coincide.  This  does  not  necessarily  eliminate  the 
possibility  of  more  than  two  “triangular"  elements  since  the  adjacent 
sides  of  a  quadrilateral  may  be  colinear  as  shown  in  figure  (6). 

(4)  The  normal  vectors  between  two  adjacent  quadrilaterals  in  a 
group  of  four  must  be  less  than  90  degrees  and  preferably  less  than  45 
degrees.  If  a  sharp  edge  is  required,  it  should  be  a  concavecornerwith 
respect  to  the  flow  field,  and  the  input  should  be  arranged  so  that  the 
edge  is  along  an  outside  boundary  of  the  groups  of  four,  and  not  through 
the  center. 

When  making  small  changes  to  the  original  geometry,  the  number  of 
elements  used  in  a  new  section  must  be  the  same  as  the  number  used  in 
the  original  section  unless  the  part  being  changed  is  at  the  end  of  the 
input  data.  Section  configurations  may  be  selected  by  natural  divisions, 
as  a  matter  of  convenience  to  more  easily  handle  large  numbers  of  points, 
or  as  a  tool  to  take  advantage  of  symmetry. 

In  setting  up  input  data  to  use  planes  of  symmetry,  it  is  important 
to  note  that  the  XYZPF  Program  has  certain  restrictions  on  the  choice  of 


symmetry  planes.  The  user  only  has  the  option  to  select  the  number  of 
symmetry  planes.  The  planes  which  will  be  used  as  symmetry  planes  are 
preselected  by  the  program  to  optimize  the  calculation  procedure. 
Therefore,  knowing  this,  the  preparation  of  input  data  must  consider  the 
following  restrictions  imposed  by  the  program  (Dean  and  Dawson  1972). 

If  only  one  plane  of  symmetry  is  used,  the  plane  of  symmetry  is  the 
y  =  0  plane  of  the  global  coordinate  system.  As  such,  all  the  y 
coordinates  of  the  input  points  must  be  of  the  same  sign,  i.  e.,  all 
positive  or  all  negative.  If  the  body  is  closed  and  intersects  the  plane  of 
symmetry,  the  points  touching  the  plane,  i.  e.,  corresponding  to  y  =  0, 
must  also  be  entered  with  the  input  points. 

If  two  planes  of  symmetry  are  used,  the  planes  of  symmetry  are  the 
y  =  0  plane  and  the  z  =  0  plane  in  the  global  coordinate  system.  Again,  the 
y  coordinates  of  all  input  points  must  have  the  same  sign,  positive  or 
negative,  and  the  z  coordinates  of  all  points  must  be  of  the  same  sign, 
positive  or  negative  without  regard  to  the  sign  of  y.  If  the  body  surface 
intersects  one  or  both  of  the  planes  of  symmetry,  the  points  which  lie  in 
the  plane,  i.  e..  those  corresponding  to  y  =  0  or  z  =  0,  must  also  be  entered 
with  the  input  points. 

If  three  planes  of  symmetry  are  used,  clearly  the  planes  are  the 
reference  planes  of  the  global  coordinate  system.  As  with  the  previous 
cases,  all  the  x  coordinates  of  the  input  points  must  be  of  the  same  sign, 
and  similarly  for  the  u  and  z  coordinates.  If  any  part  of  the  body 
intersects  any  of  the  planes  of  symmetry,  the  points  which  lie  in  that 


plane,  i.  e.,  x  =  0,  y  =  0  or  z  =  0,  must  also  be  entered  with  the  input 
points. 


5.2  SOURCE  PANEL  GEOMETRY 


Ftgore  7.  The  outer  normal  to  the  quadrilateral  element 


With  the  surface  points  identified  by  the  location  numbers,  m  and  n, 
and  arranged  in  accordancewith  program  requirements,  calculation  of 
various  aspects  of  the  source  panel  geometry  and  formation  of  the  plane 
quadrilateral  element  is  the  next  step  in  the  numerical  integration 
process.  Formation  of  all  of  the  planar  elements  is  identical,  so  the 
following  discussion  of  source  panel  geometry  will  deal  with  only  one 
characteristic  element.  The  four  cornerpoints  forming  the  basic 
quadrilateral  are  numbered  in  a  clockwise  direction  from  1  to  4  as  shown 
in  figure  (7).  It  does  not  matter  which  corner  point  is  identified  with  the 
number  1  subscript,  but  the  remaining  points  must  be  numbered 
consecutively  in  a  clockwise  direction  when  observed  from  the  flow  field 
in  order  to  ensure  an  outward  directed  normal  vector.  These  subscripts 
will  be  used  to  identify  the  cornerpoints  for  the  remainder  of  this 
discussion.  For  this  example,  the  points  will  be  identified  as  follows: 


Position  Numbers 
m,  n 
m+1,  n 
m+1,  n+1 
m,  n+1 


Global  Coordinates 
X| .  Y,  .  Z, 

X2  ,  V2  ,  Z2 
X3  ,  V3  ,  Z3 
X4  ,  V4  ,  Z4 


In  forming  the  plane  quadrilateral  elements,  the  corner  points, 
which  are  generally  not  coplanar,  are  used  to  form  the  local  coordinate 
system,  relative  to  the  element.  Recalling  that  the  crossproduct  of  two 
vectors  yields  a  vector  solution  which  is  perpendicular  to  both  of  the 
original  vectors,  the  normal  to  the  element  may  be  obtained  from  the 
cross  product  of  the  diagonals  of  the  element. 


N  =  T24  x  Tj3 


(21) 


where  T13  is  the  vector  from  point  1  to  point  3,  and  T24  is  the  vector 
from  corner  point  2  to  point  A.  The  unit  normal  is  then: 

T24  x  t13 


This  unit  normal  now  represents  the  first  of  the  three  local 


coordinate  directions,  this  one  in  the  C  direction.  The  side  of  the 
quadrilateral  from  point  2  to  point  3  is  then  used  to  obtain  the  second 
coordinate  vector. 


Ftfare  8.  The  second  local  ooordinate  vector 


N2  =  n  *  T23 


(23) 


and  the  unit  vector 


n2  = 


(24) 


Similarly,  the  third  local  coordinate  vector  is  obtained  from  the 
cross  product  of  n2  and  n,  the  result  of  which  is  a  unit  vector. 


n,z 


Fifor*  9.  The  third  tocal  coordinate  vector 


n3  -  n2  *  n 


(25) 


The  unit  vectors  n3,  n2,  and  n  form  an  orthonormal  basis  and  define 
the  local  coordinate  system  for  the  element  in  the  £,  -r\,  and  C  directions 


respectively.  Other  methods  of  obtaining  an  orthonormal  basis  could  be 
used  just  as  well,  and  would  make  no  difference  to  the  remaining 
computations.  The  origin  of  the  local  coordinate  system  would  most 
correctly  be  located  at  the  "null  point."  the  point  at  which  the  velocity 
potential  has  no  contribution  to  the  tangential  velocity  component  on  the 
source  element.  The  null  point  is  the  point  in  each  quadrilateral  element 
where  the  normal  velocity  boundary  condition  is  applied.  However,  with 
the  exception  of  long,  thin  quadrilaterals,  the  physical  difference  between 
the  null  point  and  the  centroid  of  the  quadrilateral  is  not  significant.  The 
XYZPF  Program  will  print  a  warning  in  the  output  when  a  quadrilateral  is 
long  and  thin  enough  to  jeopardize  the  accuracy  of  the  approximation  in 
that  region.  By  locating  the  origin  of  the  local  coordinate  system  at  the 
centroid,  rather  than  at  the  null  point,  the  difficult  process  of  locating 
the  null  point  for  each  element  can  be  eliminated,  later  calculations  of 
the  multipole  expansion  can  be  simplified,  and  the  boundary  conditions 
can  be  applied  at  the  centroid  without  contributing  significant  error  to 
the  approximation  (Hess  and  Smith  1966).  Therefore,  the  origin  for  each 
local  coordinate  system  is  located  at  the  centroid  for  the  respective 
element. 


5.3  LOCATING  THE  CENTROID 


The  centroid  of  the  element  may  be  calculated  by  first  dividing  the 
area  of  the  quadrilateral  into  two  triangular  areas,  the  triangles  being 
separated  by  the  line  from  point  2  to  point  4.  The  area  A^f  the  triangle 
defined  by  corner  points  2,  3,  and  4  is 

A,  =  -i-  |t24  »  T23|  (26) 

Similarly,  the  area  A2of  the  triangle  defined  by  corner  points  1,  2,  and  4 
is 

A2  =  y  |Ti2  *  Th  |  (27) 


In  the  global  coordinate  system,  the  X  component  of  the  centroid  is 
given  by 

77  A1  (28) 

A  i  +  «2 

where  X,  and  X2  are  the  averages  of  the  X  components  of  the  corner 
points  of  each  triangle.  Substituting  the  values  for  Xj  and  X2: 


W  "■ 


'-mP-rS'Jr71F^P?  *~.  r:  ^7  •V* 


vn.-'v 


x  = 


■j'Ai(X2  +  X3  +  X4)  4  A  2  ( X  ]  +  X  2  +  X  4  ) 
A,  +  A2 


3 


J_ 
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(  Aj  +  A2)X2  ♦  (  Ai  ♦  A2)X.4  4  A  ^  X3  ♦  A2Xi 
A 1  +  A  2 

A  ,  y,  +  A'-  v, 

X2  +  X4  +  ,-'a  ‘,'1 


Similarly 


Y  - 


L 


t-2  +  I4  + 


A,  ♦ 

a2  J 

.  A ,  V3 

♦  a2y, 

At 

+  A2  . 

A  1Z3 

+ 

+  A  2  Zi 

M 


1  +  A2 


(29) 


(31) 


5.4  COORDINATE  TRANSFORMATION 

Now  that  the  local  coordinate  system  is  formed  and  properly 
located  at  the  centroid  of  the  element,  the  global  coordinates  of  the 
corner  points  (X,  Y,  Z)  are  transformed  to  local  coordinates^,  r\,  0 
through  the  components  of  the  reference  vectors  of  the  local  coordinate 
system  as  follows: 


’  n?x  n?y  nlz 

X  -  X 

"  * 
s 

"2X  n2y  n2; 

V  -  V 

nx  Hg  n  j 

.  2  -  Z  . 

.  c . 

(32) 


The  corner  points  are  projected  into  the  plane  of  the  quadrilateral 
element  by  setting  the  £  components  to  zero.  The  original  diagonal 
vectors,  Ti3  and  T2«,  will  be  on  opposite  sides  of  the  resulting  plane.  The 
plane  quadrilateral  element  is  now  completely  defined.  The  program  will 
sweep  through  all  of  the  input  elements  using  the  assigned  location 
numbers,  and  repeat  this  process  for  each  element. 


Fiyare  1 1 .  Forming  the  plane  quadrilateral  element 


5.5  MOMENTS  OF  INERTIA 

The  calculation  of  the  moments  of  inertia  for  each  element  are 
performed  for  use  in  the  computation  of  the  velocity  coefficients  using 
the  quadrupole  method.  Any  calculus  text  will  give  the  moments  of 
inertia  of  a  planar  section  with  a  constant  unit  density  about  the  origin 
to  be: 


I XV  =  ’ 


B  £:'  d  t  d  T) 


(33) 


For  the  triangular  region  defined  by  the  corner  points  2,  3,  and  A, 


--  -Jy  Ui2*i3>2*(S3*i4>2*<i4*i2>2l  (36) 

I  yy  =  “  l  (  Tl2  +  )2  +  (  T| 3  ♦  ^  )2  ♦  (  ^4  +  T(2  )2  ]  (3?) 

A 

Ixy  =  77  U  42  +  £3  )(  T[2  *  %  )  +  (  £3  +  tf)(  H3  +  )  (38) 

+  (  £4  +  ^2  )(  ^4  +  ^2  )1 


Similar  equations  can  be  generated  for  the  triangular  region  defined  by 
the  corner  points  1,  2,  and  4.  The  moment  of  inertia  for  the  entire 
quadrilateral  is  the  sum  of  the  corresponding  expressions  for  each  of  the 
triangles.  The  resulting  equations  are; 

Ixx  =  I  (  i2  ♦  i3  )2  +  (  i3  *  i4)2  ♦  (  i4  ♦  i2  )2  1  (39) 

T  Hi,  *i2)2*U2*i4)2*(i4*i,  )21 

=  7  1  ( 112  *  1'3  )Z  * ( 112  *  T|4  )2  * ( 'r'i  *  ^  >2 1  (40) 

*  77  I  ( H,  *  -n.2  >2  *  (  t.2  *  ^4  >2  *  <  114  *  Tti  )Z  1 

A 

Ixy  =  7J  l(^2  +  ^)(^2  +  Tl3)  +  (^^.4)(^3  +  Tl4)  (41) 

+  (  £4  +  ^2  )(  Tt4  +  Tl2  )  ] 

+  IT  I(  ^  +  ^2)(  711  +  112  '*  +  (  ^2  +  ^4>(  *n2  +  ^4* 

+  (^4  +  ^1  K  TU  +  Tl,  )1 


>  i 


6.0  THE  MATRIX  OF  INFLUENCE  COEFFICIENTS 


With  the  quadrilaterals  completely  formed,  the  next  step  is  to 
calculate  the  velocities  induced  by  the  elements  at  the  centroids  of  all 
the  other  elements.  The  total  number  of  elements  forming  the  surface 
will  be  represented  by  N.  Let  the  source  element  be  the  (j)th  element,  and 
the  element  for  which  the  velocity  components  are  to  be  calculated  at  the 
centroid  is  the  (i)th  element.  It  does  not  matter  how  the  (i)th  elements 
are  arranged  in  relation  to  each  other  as  the  sequence  progresses. 
However,  the  sequence  must  be  consistent  as  the  calculations  proceed 
from  one  source  element  to  another.  This  program  sweeps  through  the 
(i)th  elements  in  the  order  of  their  location  numbers,  m  and  n.  For  each 
consecutiven  line,  the  elements  are  swept  in  order  of  increasing  m 
numbers. 

The  result  of  the  induced  velocity  calculations  for  the  elements 
with  unit  source  densities  is  an  N  by  N  square  matrix  of  the  values  of 
induced  velocities  at  each  element,  known  also  as  the  "matrix  of 
influence  coefficients."  The  XYZ  potential  flow  program  calculates  the 
coefficients  column  by  column,  while  the  Hess  program  calculates  them 
row  by  row.  The  advantage  of  one  over  the  other  depends  on  the  method 
of  later  solving  the  matrix  for  the  source  densities.  In  calculating  the 
influence  coefficients,  twenty-five  quantities  which  describe  the 
geometry  of  the  source  element  are  required  to  adequately  calculate  the 
induced  velocity  at  the  centroid  of  the  (i)th  element.  These  quantities 
include  the  coordinates  of  the  centroid  in  the  global  coordinate  system, 
the  elements  of  the  coordinate  transformation  matrix,  the  local 
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coordinates  of  the  corner  points,  the  maximum  diagonal,  the  area,  and  the 
second  moments  of  the  quadrilateral  element.  Additionally,  the  Hess 
program  uses  the  coordinates  of  the  null  point,  making  a  total  of 
twenty-eight  quantities  for  that  method  (Hess  and  Smith  1962). 

When  calculating  row  by  row,  the  first  (i)th  element  is  selected, 
containing  the  "null"  point,  and  the  influence  coefficients  are  computed 
for  all  of  the  (j)th  elements  in  sequence  before  proceeding  to  the  (i+1)th 
element.  This  procedure  requires  the  twenty-five  quantities  for  each 
(j)th  element  to  be  available  for  calculation  of  the  influence  coefficient. 
Because  each  of  the  N  (j)th  elements  is  used  N  times  with  this  procedure, 
calculating  the  geometric  quantities  or  retrieving  the  values  from  low 
speed  memory  would  be  very  time  consuming,  since  the  calculations  or 
memory  access  would  need  to  be  performed  N2  times.  Therefore,  it  is 
more  practical  to  have  the  values  available  in  high  speed  memory, 
although  large  matrices  may  exceed  the  storage  capacity  of  high  speed 
memory,  imposing  a  limit  on  the  number  of  elements  which  can  be  used. 
The  advantage  to  the  row-by-row  calculation  is  that  solution  of  the 
resulting  matrix  by  the  Gauss-Seidel  reduction  method  does  not  require 
transposing  the  matrix,  which  would  be  another  time  consuming  process 
(Hess  and  Smith  1962). 

Another  alternative  is  calculation  of  the  influence  coefficients 
column  by  column.  This  method  calculates  the  influence  coefficients  by 
sweeping  all  the  (i)th  elements  for  each  (j)th  element  before  proceeding 
to  the  (j+1)th  element.  Therefore,  the  twenty-five  geometric  quantities 
are  retrieved  from  low  speed  storage  only  once  for  each  (j)th  element, 


for  a  total  of  N  times.  This  procedure  is  not  limited  by  the  capacity  of 
high  speed  memory,  and  calculation  of  the  coefficient  matrix  is 
approximately  10&  faster  than  the  row-by-row  method  (Dawson  and  Dean 
1972).  This  is  the  calculation  method  used  by  the  XYZ  Potential  Flow 
Program. 

An  influence  coefficient  represents  the  combined  effects  on  one 
element  of  the  velocity  potentials  of  all  the  other  elements  comprising 
the  body  surface.  For  the  quadrilateral  element  with  a  unit  source  density 
in  the  xy-plane,  from  equation  (13),  the  potential  at  point  P  (x,  y,  z)due 
to  the  element  is 


The  integrand,  l/r,  can  be  expanded  in  a  series  about  the  origin  in 
powers  of  £  and  t\.  Each  term  of  the  series  will  contain  the  product  of 
some  powers  of  £  and  t|  with  a  corresponding  derivative  of  l/r0,  where  r0 
is  the  distance  of  the  field  point  P  from  the  quadrilateral  origin. 


Then  the  series  expansion  through  the  second  order  term  is  (Hess  and 
Smith  1962): 
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<P  -  Aw  -  (Mxwx +  ^ywy) +  ^(I  xxwxx  +  ^  *xywxy +  *yywyy^ +  •  •  •  (^) 

The  subscripts,  x  and  y.  used  with  w  represent  the  respective  partial 
derivatives.  This  series  represents  the  multipole  expansion  of  the 
velocity  potential,  since  each  term  can  be  interpreted  as  a  point 
singularity  of  a  particular  order.  The  first  term  is  the  potential  at  point 
P  due  to  a  point  source  of  strength  A  located  at  the  origin.  The  second 
term  is  the  sum  of  two  dipoles  of  strengths  Mx  and  My  located  at  the 

origin,  oriented  along  the  x  and  y  axis  respectively.  The  choice  of  the 
centroid  of  the  quadrilateral  as  the  origin  of  the  local  coordinate  system 
causes  the  first  moments,  Mx  and  My,  to  be  zero.  Therefore,  the  dipole 
terms  disappear,  and  are  not  dealt  with  anywhere  in  the  program.  The 
third  term  is  the  sum  of  three  quadrupoles  of  strengths  Ixx,  lXy,  and  lyy 

located  at  the  origin.  Kellogg  (1929)  shows  that  this  secondorder 
approximation  is  absolutely  and  uniformly  convergent,  and  Hess  and  Smith 
(1962)  show  that  convergenceis  rapid  enough  with  an  increase  in  r0  that 
certain  further  approximations  may  be  made  without  significant  error  at 
large  distances  r0  from  the  source  quadrilateral. 

Hess  and  Smith  (1962)  presented  a  comparison  of  velocities 
calculated  using  the  exact  formulas,  a  simple  point  source,  and  a  second 
order  approximation.  The  comparisons  were  based  on  the  ratio  of  the 
distance  r0,  between  the  centroid  of  the  source  quadrilateral  and  the  field 
point  P,  to  the  length  of  the  maximum  dimension  t,  of  the  source 
quadrilateral,  typically  the  maximum  diagonal.  The  non-dimensional  ratio 
is  then  r0/t.  The  results  show  that  sufficient  accuracy  is  maintained 
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while  using  a  simple  point  source  at  ratios  of  (r0/t)  >  4,  using  the  second 
order  source  and  quadrupole  solution  for  the  range  2.45  <  (r0/t)  <  4,  and 
using  the  exact  solution  for  ratios  of  (r0/t)  <  2.45.  In  any  case,  the  error 
goes  to  infinity  as  the  field  point  approaches  the  edge  of  the  quadrilateral 
where  calculations  indicate  an  infinite  velocity.  TheXYZ  Potential  Flow 
Program  uses  a  monopole  source  for  (r0/t)  >  4,  the  source  -  quadrupole 
formulae  for  2  <  (r0/t)  <  4,  and  the  exact  formulae  for  (r0/t)  <  2.  Hess 
and  Smith  (1962)  reported  a  maximum  error  of  0.001  in  approximating  any 
velocity  component  using  the  above  criteria. 


From  equations  (2)  and  (42).  the  components  of  the  velocity  at  the 
field  point  p(x.  y.  z)due  to  the  source  quadrilateral  are; 


_  ft  (X  -  j)  d£  dn 

JJ  r3 


f  (y  -  T[)  d£  dTj 

1  r  3 


89  Cf  z  d£  d'n 

r  ft  r 3 


Equations  (44).  (45)  and  (46)  are  evaluated  by  expressing  each  of 
the  integrals  as  the  sum  of  four  terms,  each  term  representing  the  effect 
of  one  side  of  the  quadrilateral  (Hess  and  Smith  1962).  This  method  can 
also  be  generalized  for  polygonal  elements  with  any  number  of  sides.  The 
potential  function  for  each  side  of  the  quadrilateral  is  the  combined 
potentials  of  semi-infinite  strips  whose  boundaries  are  the  side  of  the 
quadrilateral  and  two  semi-infinite  lines  parallel  to  either  the  x  or  y 
axis.  When  observed  from  the  domain,  and  the  sides  are  traversed  in  a 
clockwise  direction,  the  source  strip  on  the  right  will  have  a  source 
density  of  a  =  +1/2  and  the  sourcestrip  on  the  left  will  have  a  source 
density  of  c  =  -1/2  as  shown  in  figure  (12).  When  the  sides  are 
recombined  to  form  the  quadrilateral,  the  source  densities  outside  the 
quadrilateral  cancel  each  other,  and  the  source  densities  within  the 
quadrilateral  combine  to  form  a  source  density  of  O  =  +1.  This  will  be 
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true  for  a  planar  element  with  any  number  of  sides  and  in  any  relative 
orientation  within  the  plane. 


7.1  THE  Y  VELOCITY  COMPONENT 

From  equation  (45),  the  velocity  component  Vy  is  found  by  summing 

the  four  terms  representing  the  contributions  of  the  sides  of  the 
quadrilateral.  For  the  side  from  point  (£,,  ti,)  to  point  (£2,  t}2),  the 
contribution  is  expressed  as  the  integral  over  the  area  of  the 
semi-infinite  strips  with  the  source  densities  of  o  =  +1/2  and  o  =  -1/2 
rather  than  the  unit  source  density  of  equation  (45). 
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Integrating  with  respect  to  -q: 


V  J2  2 
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The  terms  evaluated  at  tj  =  +«>  and  t\  =  -»  cancel,  and  the  terms  evaluated 
at  t\  =  t|12  add  to  obtain  the  following  expression-. 

v'2  =  2  J  t,  Kx-i)2*(a-T()2*22]H 

vy,£  =  J  „  —  (48) 

J  <i 

Equation  (48)  is  changed  to  a  function  of  arclength  s  by  the  relation 


\!  ( s2 ~  hr  +  ^2- V 


where  d12  is  the  length  of  the  side  of  the  quadrilateral  fromCS,.  t\,)  to 
(£2,  t\2)  as  shown  in  figure  (13). 


Fifrt  1 3.  The  potential  due  to  a  finite  tine  source  (Hess  &.  Smith  1  %2) 


Substituting  equation  (49)  into  equation  (48) 
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y 


dl2  Jo  r 


(50) 


From  figure  (13),  it  can  be  seen  that,  in  terms  of  arclength  s,  the  distance 
r  from  point  P  to  any  point  on  the  line  from  point  1  to  point  2  is  given  by 


r 


=  'n  fz  +  (si  -  i  'r 


Substituting  into  equation  (50)  yields 


Evaluating  the  integral 
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(52) 


The  quantities  r,f  r2.  S\,  and  s2  used  in  equation  (52)  are  as  shown  in 
figure  (13).  Equation  (52)  is  singular  when  r,  =  s,,  which  occurs  when  the 
field  point  P  is  located  anywhere  along  the  line  defined  by  the  side  of  the 
quadrilateral.  This  singularity  may  be  removed  by  using  the  law  of 
cosines  (Hess  and  Smith  1962). 


From  equation  (52) 


y,£  ~  d12  ’~3  f  1  ( 1  -  COS  £<  ) 

where  and  fi2  are  interior  angles  shown  in  figure  (13).  Applying  the 
law  of  cosines  to  figure  (13) 
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From  equations  (54)  and  (55): 
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Substituting  equation  (56)  into  equation  (53)  yields  the  final  form  of  the 
exact  equation  of  the  y  component  of  the  velocity  induced  by  the  side  of 
the  quadrilateral  from  point  1  to  point  2; 


Equation  (57)  is  applied  to  the  remaining  sides  of  the  quadrilateral 
simply  by  substituting  the  appropriate  point  numbers  for  the  corner 
points  of  each  side.  The  total  contribution  of  the  quadrilateral  to  the  y 
component  of  the  velocity  is  the  sum  of  the  four  terms  representing  the 
contributions  of  each  of  the  sides.  The  y  component  of  the  velocity  at 
the  field  point  P  is  now  given  by: 


7.2  THE  X  VELOCITY  COMPONENT 


A  similar  derivation  process  is  used  to  produce  the  equation  for  the  x 
component  of  the  velocity  induced  by  the  side  of  the  quadrilateral  from 
point  1  to  point  2.  The  semi-infinite  source  strips  are  constructed 
parallel  to  the  x  axis,  and  the  order  of  integration  is  reversed. 
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The  x  component  of  the  velocity  at  the  field  point  P  due  to  the 
quadrilateral  is  given  by: 


d12 

^4  -Tl* 


Inn 

r,  +  r2  +  d|2 

%  -  ^2 

1  fifl 

r2  +  r3  +  d23 

log 

r,  +  r2  -  d|2 

d23 

1  U 

r2  +  r3  ~  d23 

Hs  +  r4  +  d34 

Tii  -  n4 

1  nn 

r4  +  n  +  d4i 

lug  • 

r3+  r4  -  d34 

d4i 

lug 

r4  +  n  _  d41 

7.3  THE  Z  VELOCITY  COMPONENT 


The  z  component  of  the  velocity  at  the  field  point  P  due  to  the 
quadrilateral  is  obtained  in  a  similar  fashion,  using  semi-infinite  source 
strips,  this  time  parallel  to  the  y  axis.  From  equation  (46),  the 
fundamental  velocity  potential  of  the  semi  infinite  source  strips  is 
integrated  in  a  manner  similar  to  that  used  to  obtain  equation  (47),  and 
the  z  component  of  the  velocity  due  to  the  side  from  ($,,  7|t)  to  (£2,  t[2)  is 
given  by 

V2r.  =  T  d$  [i  f  t|2"  ^  P  ]  — — —  (60) 

J  E;,  L "  J- co  ‘  J  71,2  J  r3 

,,  =  z.  f*2  [  C'2-  ri  dn _ 

Z'2  2  7  E,  i-J-w  A,J[(X  -  ^)2  +  (y  "  Tl^  +  Z8]54 


Performing  the  integration  with  respect  to  t\,  the  integral 


fits  the  integral  form 


where 


[C2+F2]n  2Cnn-1)  [[C2  +  F*-]' 

i 

C2  =  (x  -  O2  ♦  z2 
F  =  (y  -  t\) 
dF  =  -dn 
n  =  3/2 


[C2+  F2r] 


Then,  from  equation  (60) 


__z  r _ (y  -  Tj) _ 

2  J  x,  u'""  i  Kx  -  o2  *  2 2  ][  (x  -  Oz  +  (y  -  ti)£  +  z2 nia, 
(y  -  'a) 

Kx-£)£  +  2£][(x-£)£  +  (y-Ti)£  +  2£]>2  _co 

_ (y  -  T|) _  (& 

Kx  -  £)£  +  2 £ ][  (x  -  02  +  (y  -  ti)£  +  22]^ 

_ (y-T\) _  j  ^ 

[(x  -  t.)2  +  z 2  ][  (x  -  £)£  ♦  (u  -  i\)2  +  z£  J'‘  L  / 


Again,  the  terms  evaluated  at  +°°  and  cancel  and  the  terms  evaluated  at 
t\12  add  to  obtain  the  following  expression: 


(y-WdS 


Kx  -  $)2  +  z2H(x  -  O  +(y  -  ti,2) 


>: 

**. 

? 

•s' 

V 


u 


Without  a  convenient  substitution  with  which  to  integrate  equation  (62), 
the  integration  is  performed  directly.  Recognizing  that  along  the  line 
defined  by  the  side  of  the  quadrilateral  from  (£,.  *q,)  to  (£2,  t\2).  t\\2  may 
be  expressed  as  a  function  of  £: 


Tli2  =  m,2  £  +  b,2 

where  the  s’ope  of  the  side,  m12  is  given  by 


m,2 


S2  -  Si 


(63) 


(64) 


and  b12  may  be  determined  knowing  that  t}|2  =  tj,  when  S  =  Si- 


!I2 


Substituting  equation  (63)  into  equation  (62)  yields 


(.65) 


■ 


r  _ (y  _ 

J  Kx  -  S)2  +  z2)((x  -  S)2  +  (y  -m!2S-b,2)2  +  z23^ 


(66) 


Define  the  quantities 


qi2  =  y  -  b]2  -  m12x  (67) 

u  =  x  -  S  (68) 

Then  du  =  -d£  (69) 

y  -  b12  -  m12S  =  m12u  ♦  q12  (70) 


By  a  change  of  variable,  equation  (66)  is  now  expressed  as  a  function  of  u 


which  fits  the  form  of 

r _ (lu  +  n)  du _ 

J  (Au2  +  2Bu  +  C)4  (au2  +  2bu  +  c)‘ 

where 


L  =  m12 

H  =  qi2 

A  =  1 

a  =  (m122  +  0 

B  =  0 

b  =  m12q12 

C  =  z2 

c=q122  +  z2 

From  Hardy  (1944),  this  integral  form  may  be  integrated  by  the 
substitution 

...  I't  +  q.  (74) 

t  +  1 

where  p  and  u  satisfy 

ajju  +  b(p  +  u)  +  c  =  0  (75) 

Ajiu  +  B(y  +  u)  +  C  =  0  (76) 

and  are  the  roots  of  the  equation 


(aB  -  bAK2  -  (cA-  aC)  *  ♦  (bC  -  cB)=  0 


(77) 


Substituting  the  appropriate  values  into  equation  (75),  the  roots  of  the 
quadratic  equation  are 


M  =~  — 

m12 

(?8) 

m,2z* 
v  =  — - — 

(79) 

it  can  be  verified  that  these  values  satisfy  equations  (75)  and  (76). 

Substituting  equations  (78)  and  (79)  into  equation  (74) 

m12z2  q12t 

Qi2  ‘  m12 

t  ♦  1 

(80) 

r  m,2z2  -I 

du  =  '  m’2  dt 

(t  +  I)2 

(61) 

By  substitution  and  a  change  of  variable,  equation  (72)  becomes  a  function 

of  the  parameter  t.  After  simplification,  the  integral  now  fits  the  form 

of 

r  ...  ..  0t 

K  J  (c<t2+  3)4(#t2  +  5)’ 

(82) 

where 

K  =  -(m,26q12z*  ♦  2m124q123z2  +  qi25m122) 

«  =  0i24  +  rn122q)22z2 

50 


3  =  m124z4  ♦  m122q122z2 
*  =  qi24  ♦  m122q122z2 

8  =  m126z4  +  m124z4  +  2m124q122z2  +  m122q124  + 

m122q122z2 

Equation  (82)  can  be  rationalized  by  the  substitution 

t 

(  O  T 

•sftfFTT 

KOo) 

from  which  it  can  be  shown  that 

v!£ 

1  -  v2  '6 

(R4) 

8  1V2 

dt=  {.  (1  -v*tf)9  J  0v 

(85) 

Substituting  equations  (84)  and  (85)  into  equation  (82)  and  simplifying 

yields  the  integral  in  terms  of  the  parameter  v: 

J  (Kt2  +  £>J(tft2+  8) 

which  fits  the  form 

r  t>v 

J  a2  +  b4 v£ 


$  +  (c<  6  -  3  t5 )  v 


1  ,  -i  Dv 

ih  tan  T 


where  a2  =  3 

t>2  =  (o<S  -  fit) 


Performing  the  integration 


4£(<xS-£b) 


£  J 


From  equations  (80),  (83),  and  the  expressions  fora,  fi,  5,  and  2f  from 
equation  (82),  and  after  a  considerable  amount  of  algebraic  manipulation 
and  simplification,  equation  (88)  becomes 


4  6  -  $  b  ) 


tan1  v 


tan"1 


'  miaza-  gig  u 
z  Jza  +  ua  +  (ml2u  +  q«)2‘ 


From  equations  (63)  and  (67),  equation  (89)  becomes 


2  +  u c'  +  (miau  +  412, 


C  m,2(u2  +  z2) 


Finally,  applying  these  results  to  equation  (71) 


(x  -  Ci  )2+(y  -  T(12)‘ +  z£ 


2  ( (x  -  6  )2  +  z2)  -  (y  -  ti12)(x  -6)1 


2  4(x  -  6  )e+  (y  -  T|12r  +  z2 

Recall  that  when  x  =  6,  y  =  'Hi  and  when  x  =  £2.  y  =  ti2.  Then,  for  the  sake 
of  a  more  compact  equation,  define  the  following  quantities-, 
e,  =  (x  -  6)2  +  z2  e2  =  (x  -  £2)2  +  z2 

hi  =  (y  -  tiOCk  -  6)  h2  =  (y  -  *q2)(x  -  6) 

The  quantities  r,  and  r2  are  as  shown  in  figure  (12),  where 


r1  =  4(x-6)a  +  (y-'n1)a  +  za 


r2  =  4(x  -  6)£+  (y  -  ti£)£  +  z£ 


Substituting  these  quantities  into  equation  (91)  yields  the  exact  z 
component  of  velocity  due  to  the  side  from  point  (6/Hi)  to  (£2,t[2)  in  the 
form  used  by  the  XYZ  Potential  Flow  Program: 


[  m12ei  -  hi  1 


r  m,2  e2  -  h2 1 


I 


» 


L 


i 


8.1  QUADRUPOLE  METHOD 

As  previously  mentioned,  as  the  ratio  of  (r0/t)  exceeds  the  value  of 
2,  then  certain  approximations  may  be  made  which  greatly  reduce  the 
calculation  effort  otherwise  required  by  the  exact  method.  In  the  range 
of  2  <  (r0/t)  <  4,  the  XYZ  Potential  Flow  Program  uses  the  second  order 
approximation  of  the  potential  described  by  equation  (43).  With  the 
origin  at  the  centroid  of  the  quadrilateral,  the  first  moments  are  zero, 
and  the  second  order  approximation  is 

<p  =  Aw  *  (l/2)(!  xxwxx  +  2,xywxy  +  ^yywytP  (94) 


where  the  first  term  is  a  point  source  of  strength  A,  the  second  term  is 
composedof  three  quadrupoles  of  strengths  lxx,  IXy .  and  lyy  located  at 

the  local  origin,  and  the  subscripts  onw  indicate  the  partial  derivatives 
of  w  with  respect  to  those  variables  as  before.  The  quantity  A  is  the  area 
of  the  element,  and  the  terms  Ixx,  lxy,  and  lyy  are  the  respective  moments 

of  inertia  of  the  source  element  given  by  equations  (39),  (40),  and  (41). 

To  obtain  the  velocity  components  at  the  field  point,  equation  (94)  is 
differentiated  with  respect  to  the  coordinate  directions  giving: 


V  X-  =-[AWx+  2  ^xx  WXXX+  Ixy  WXXy  +  2  lyy  WXyy  J  (95) 

Vy  =  ~  |^-  =  -[  AWg  +  4  IxxWxxy+  IxyWXyy+  ^IyyWyyy]  (96) 

Vz=”dF  =  -  [AWZ  +  \  IXXWXX3+  IxyWXyZ+  -jlyyWyyz]  (9?) 

Recallinq  that  w  =  ,  -  ■  =r  =  7- 

-I  x2  +  y2  +  2 2  r° 


the  derivatives  of  w,  as  expressed  by  Hess  and  Smith  (1962)  and  as  used  in 
the  XYZPF  program,  are 


VVj;  = 

-x  r0-3 

Wy  = 

-y  ro~3 

> 

YVt  — 

-zr0-3 

W 

XXX 

=  3x(3p  + 

1 0x2) 

VYxxy 

=  3y  p  ty' 

W.. 

xyy 

II 

X 

-O 

-1 

VVyyy 

=  3y(3q  + 

10y2) 

WXXZ 

=  3z  p  y" 

Wxyz 

=  - 1 5  x  y  z 

y7 

wyyz 

=  3z  q  y- 

7 

.  where 


( 1 00) 


p  =  y2  +  z2  -  4x2 
q  =  x2  ♦  z2  -  4y2 
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8.2  MONOPOLE  METHOD 


When  the  ratio  of  (r0/t)  is  greater  than  4,  then  the  quadrilateral 
may  be  approximated  by  a  simple  source  corresponding  to  the  first  term 
of  equation  (43).  Then  the  velocity  components  at  the  field  point  due  to 
the  quadrilateral  are  given  by 


---AW* 


(101) 


v*  =  -  %  =  -  A"'n 
Vj=  ~w  -AW, 


( 1 02) 


(103) 


where  the  partial  derivatives  of  w  are  those  given  in  equation  (99). 


v  W  Si 


s.s 


9.1  JACOBI’S  ITERATIVE  METHOD 


From  equation  (20),  the  matrix  equation  may  be  solved  for  the 
constant  source  density  Oj  for  each  element  which  satisfies  the  boundary 

condition  equation  (II).  Equation  (20)  suggests  the  use  of  Jacobi’s 
iterative  method  of  matrix  solution  in  the  form 


(m+1) 


Vf  *  Y.  •  '  =  1.2 . N 


(104) 


where  N  is  the  number  of  elements  composing  the  body  surface,  and  m  is 
the  number  of  iterations  completed.  A  partial  sum  of  equation  (104)  is 
computed  for  each  of  the  ith  elements  before  proceeding  to  the  next  jth 
element.  The  iteration  is  complete  when  the  summation  of  equation 
(104)  includes  all  of  the  jth  elements.  Because  the  values  of  the  source 
densities  at  all  of  the  elements  are  recomputed  before  any  of  them  are 
used  in  the  iteration,  this  method  is  also  called  the  simultaneous 
displacement  method  (Ralston  1965).  This  is  contrasted  with  the 
Gauss-Seidel  iterative  method  used  in  the  Douglas  program.  In  the 
Gauss-5eide!  method,  as  each  new  Oj  is  computed,  it  is  used  immediately 

in  the  iteration  process  for  calculation  of  0(i+])  •  This  is  also  known  as 
the  successive  displacement  method  and  is  expressed  as 


.N  A  **  V 


■>>  W  J>Vv 


V.  *  Z  Cyor^i  C,;0«  (105) 

j=1  j=H1 

j  =  1,  2,  ....  N 

Though  the  Gauss-Seidel  iterative  method  is  faster,  the  Jacobi  iteration 
method  was  selected  for  use  in  the  XYZPF  program  in  order  to  be  able  to 
perform  the  iterations  column  by  column,  since  the  coefficient  matrix  is 
also  computed  column  by  column,  and  the  matrix  does  not  have  to  be 
transposed  for  solution. 

When  the  (m+1)th  iteration  is  complete,  the  values  of  the  source 
densities  are  compared  with  those  of  the  (m)th  iteration  and  the 
differences  summed  for  all  of  the  elements.  The  total  difference  between 
successive  iterations  is  then  compared  to  a  convergencecriteria  input  by 
the  user.  If  the  difference  is  less  than  the  convergencecriteria,  then  the 
matrix  solution  is  complete  and  the  values  of  the  source  densities  are 
stored  for  later  use  in  computing  velocities  and  pressure  coefficients.  If 
the  convergencecriteria  is  not  met,  then  the  iteration  process  is 
repeated.  After  every  five  iterations,  if  the  convergencecriteria  is  still 
not  met,  then  an  extrapolation  is  attempted  in  order  to  accelerate  the 
convergence.  The  XYZ  Potential  Flow  Program  uses  a  Richardson 
extrapolation  method,  a  numerical  procedure  which  uses  two  approximate 
results  to  obtain  a  third  approximation  which  is  closer  to  the  exact 


(m+1) 


solution  (Ralston  1965). 
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9.2  RICHARDSON  EXTRAPOLATION 


The  Richardson  extrapolation  assumes  that  the  iterative  process  is 
convergent.  For  the  iterative  solutions  Sq.  S1(  and  S2.  where  Sq  is  the 
most  recent  approximation  and  S2  the  oldest,  the  solution  is  convergent 


5q  -  Si 
Si  -  s2 


=  X  < 


(106) 


While  a  Richardson-type  extrapolation  can  take  may  forms,  the  XYZPF 
program  uses  a  procedure  developed  from  the  following  approximations 
(Dawson  and  Dean  1972).  If  there  is  only  one  dominant  eigenvalue  and  a 
sufficient  number  of  iterations  have  been  completed,  the  iterative 
solutions  may  be  approximated  by 


S0  *  Sf  +  E  Xm 
S,  *  Sf  ♦  E  Xm_1 
S2  ~  Sf  +  E  Xm"2 
Sj  %  Sf  ♦  E  Xm_i 


007) 


where  is  the  true  solution 

X  is  the  eigenvalue 

E  is  the  eigenfunction 

m  is  the  number  of  completed  iterations 

Define  the  linear  combination  which,  from  equation  (107),  may  be 
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approximated  as 

ASq  +  (1  -  A) S,  Sf  ♦  E  Xn"j  (AX  ♦  1  -  A)  (108) 

The  value  of  A  may  be  chosen  such  that 

(AX+1-A)=0  (109) 

Then,  from  equations  (108)  and  (109) 

AS0  +  (1  -A)  5,  *Sf  (110) 

where  the  expression  on  the  left  converges  to  the  exact  solution. 
From  equations  (106)  and  (109) 


5q  -  Si 
Si  -  S2 


Solving  for  A. 


(!!!) 


a  -  S2  ~  Si  _  S2  ~  Si  ( j  j?) 

Sq  -  2Si  +  S-2  ’  D 

Since  the  value  of  A  generally  changes  from  element  to  element,  a 
weighted  average  of  A  is  used  in  the  extrapolation,  where 


A  = 


^  (S2(0-Si(i))(sign  of  D(i)) 
1=1 


N 

t  w>) 

i=1 


( 1  1 3) 


Equation  (113)  is  recomputed  after  every  fifth  iteration.  If  the  difference 
between  the  new  value  and  the  old  value  is  less  than  0.02,  then  the 
solution  is  extrapolated.  From  equation  (110).  the  extrapolated  solution 


">■">  -*■  J  j  .*  -v  r.»  v  r 


When  there  are  two  dominant  eigenvalues,  then  the  iterative  solutions 
may  be  approximated  by 

Sj  «  Sf  ♦  E,  X,m-j  +  E2  X2m_i  (1,5) 

where  Sf  is  the  true  solution 

X,  and  X2  are  the  eigenvalues 

Ej  and  E2  are  the  eigenfunctions 

m  is  the  number  of  completed  iterations 

Define  the  linear  combination  which,  from  equation  (115),  may  be 
approximated  as 


B2  Sg  +  +  (1  ~  B^  -  B2)  S2 

58  5f  +  E]  X,"*"2  IB2  X,2  +  B,  X,  +  (1  -  B,  -  B2)]  (116) 

+  E2  X^2  [B2  X22  ♦  B,  X2  ♦  (1  -  B,  -  B2)l 

The  values  of  B,  and  B2  may  be  determined  for  which  the  eigenvalues  X, 
and  X2  are  roots  of  the  quadratic  equation 

B2  X2  ♦  B,  X  ♦  (1  -  B,  -  B2)  =  0  (117) 

Then,  from  equation  (116) 


B2  5g  +  B(  5^  +  (l  ~  B^  ~  B2)  52  ^  5f 


018) 


where  the  left  side  of  the  equation  (118)  convergesto  the  exact  solution. 
Using  equation  (115)  and  eliminating  terms  containing  E2: 


(S0  -  S,)  -  X2  (S,  -  S2)  =  E,  Ar2(A,  -X2)(X,  -  I) 

(s,  -  S2)  -  X2  (S2  -  S3)  =  E,  xr3(x,  -X2)(X,  -  1)  (119) 

(S2  -S3)  -  X2  (S3  -  S4)  =  E,  X,ft-4(X1  -X2)(X,  -  I) 

Solving  for  X, 

}  ,  -  (5o~  1  ~  X2O1  -  5a)  _  (Si  -  S2)  -  (5a  -  S3)  ^  1 20) 

(Si  -  Sa )  -  X2  (S2  -  S3 )  ( S 2  -  S3 )  _  X2  (S3  -  S4) 

From  equations  (117)  and  (120) 

_  (S4-  S3)(So-2S2+S4)-(S4- S2)l(Si  -  S2MS3- S4)l  ,  , 

e,  = - - - —  (121) 

g  (S-*~  52)(S4-2S3^S2)-(S4-  5s) [(Si  -  S2)-(Ss-  S4)]  (]^o) 

where  D  =  (S4  -  2S3  -  S2)(S0  -  2S2  +  S4)  -  (S,  -  S2  -  S3  +  S4)2 


The  weighted  averages  of  Bj  and  B2  are  used  for  the  extrapolation  as  done  i 

with  A  in  equation  (113).  If  the  sum  of  the  absolute  values  of  the  ! 

weighted  averages  of  Bt  and  B2  changes  by  less  than  2%,  then  the 

extrapolation  is  performed.  Then  from  equation  (118),  the  extrapolated 

solution  is  ! 

s*=  si  So +  eisi+  ( f  -  si  -§2)52  023)  : 

l 
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10.0 


With  the  influence  coefficients  and  the  source  densities 
determined,  the  calculation  of  velocities  is  a  relatively  simple  matter. 
From  equation  (9),  the  total  velocity  is  the  sum  of  the  freestream 
velocity  and  the  disturbance  velocity  due  to  the  body.  The  product  of  the 
source  densities  and  the  influence  coefficients  are  summed  for  all  of  the 
elements,  and  then  added  to  the  freestream  velocity  to  determine  the 
total  velocity  at  any  point  in  the  domain.  Velocities  on  the  surface  of  the 
body  are  calculated  at  the  null  points  only,  as  the  boundary  conditions  are 
enforced  only  at  the  null  point  of  each  element,  and  velocities  at  other 
points  in  the  element  would  produce  significant  error  due  to  the  method 
of  approximation.  The  components  of  the  velocity  at  the  centroid  of  the 
ith  element  are 


From  equation  (15),  the  velocity  induced  by  an  element  at  its  own  null 
.  point  has  a  magnitude  of  2tt  directed  along  the  outward  normal  vector  of 
the  element. 
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At  a  point  off  the  surface  of  the  body,  the  components  of  the 
velocity  are  determined  just  as  if  the  point  of  interest  was  a  null  point  of 
a  single  element.  The  total  velocity  at  the  field  point  is  the  sum  of  the 
freestream  velocity  and  the  contributions  of  each  of  the  elements  of  the 
body  surface.  The  contribution  of  each  of  the  elements  is  determined  by 
calculating  the  influence  coefficient  based  on  the  element  geometry,  and 
multiplying  the  result  by  the  source  density  for  the  element.  The  total 
velocity  at  the  field  point  may  be  expressed  as 

N 

Vp  =  Vw  +  ^  Cpq  Oq  (125) 

<1=1 


where  p  and  q  represent  the  field  point  and  the  source  element 
respectively  and  the  influence  coefficient, 


As  discussed  in  section  6.0,  the  influence  coefficient  may  be 
calculated  by  the  exact  method,  or  it  may  be  approximated  by  the 
quadrupoie  or  monopole  method  depending  on  the  ratio  of  the  distance,  r0, 
between  the  field  point  and  the  centroid  of  the  source  element  to  the 
maximum  dimension,  t,  of  the  source  element. 


The  magnitude  of  the  velocity  at  either  the  on-body  or  off-body  ; 

i 

points  is  given  by  ■ 

i 


027)  ; 

i 

I 
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The  pressure  coefficient  is  calculated  by  using  the  result  of 


I 


1 
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For  the  calculation  of  streamlines  off  the  surface  of  the  body,  a 
timestep  procedure  is  performed  by  calculating  the  velocity  at  the 
starting  point  of  the  streamline  from  equation  (125),  and  advancing  the 
streamline  onetime  increment  by  a  fourth  order Runge-Kutta  integration 
to  a  new  point  (Ralston  1965).  The  timestep  procedure  is  repeated,  thus 
creating  a  streamline  composedof  finite  segments. 


For  the  calculation  of  streamlines  on  the  surface  of  the  body,  the 
streamline  is  started  at  a  specified  point  and  quadrilateral  number.  The 
local  velocity  is  calculated  from  equation  (124),  and  the  values  of  a 
stream  function  are  computed  for  each  corner  point.  The  stream  function 
is  chosenso  that  it  has  a  value  of  zeroat  the  last  point  on  the  streamline 
in  the  quadrilateral.  The  side  of  the  quadrilateral  through  which  the 
streamline  exits  is  determined,  and  coordinates  of  the  point  on  the  side 
which  has  a  stream  function  value  of  zero  are  computed.  The  direction  of 
the  streamline  is  verified  by  comparing  it  with  the  known  direction  of 
positive  velocity.  The  next  quadrilateral  through  which  the  streamline 
passes  is  determined  by  calculating  the  proximity  of  the  new 
quadrilateral  to  the  most  recent  point  on  the  streamline.  A  circular  area 
is  computed  which  encloses  the  new  quadrilateral  with  an  additional  10& 
margin.  If  the  last  point  of  the  streamline  falls  outside  the  circle,  then 
the  quadrilateral  is  discarded  and  a  new  one  selected  until  the  streamline 
is  adjacent  to  the  new  quadrilateral. 


Ftfar*  IS.  Calculation  of  on-body  stroamHnes 


This  procedure  is  repeated  along  the  surface  of  the  body  until  all  of 
the  surface  elements  have  been  tested.  The  result  is  a  streamline 
composed  of  segments  from  one  side  of  an  element  to  another. 


12-0 


The  XYZ  Potential  Flow  program  assumes  a  constant  element 
source  pane!  as  described  in  Section  3.3.  Extensive  use  of  the  constant 
element  source  panel  method  has  shown  that  the  primary  diasadvantage  of 
the  method  is  that,  in  order  to  obtain  a  highly  accurate  solution,  a  large 
number  of  surface  elements  must  be  used  to  discretize  the  body  surface. 
The  method  has  been  applied  to  problems  of  increasingly  complex 
configurations  (Hess  1977).  By  doing  so,  the  size  of  the  coefficient 
matrix  is  increased  resulting  in  increased  computer  time  and  cost. 
Additional  cost  is  accrued  due  to  the  manhours  required  to  prepare  the 
input.  Therefore,  while  the  constant  element  methods  have  proven  to  be 
very  successful,  the  cost  has  motivated  the  development  of  higher  order 
methods. 

The  higher  order  surface  singularity  methods  discretize  the  body 
surface  with  curved  elements  having  a  variable  source  density,  as 
compared  to  the  flat  elements  of  constant  source  strength  used  in  the 
basic  method.  Hess  (1973)  showed  that  the  effect  of  a  curved  surface  and 
the  effect  of  a  variable  source  density  are  of  the  same  order  of 
magnitude.  Therefore,  the  two  effects  must  be  used  together  to  provide  a 
“consistent"  solution.  The  consistent  higher  order  panel  method  provides 
the  increased  accuracy  and  speed  desired  for  three  dimensional  Neumann 
-  problems  (Hess  1979). 

Accordingto  Hess  (1979),  the  evolution  of  the  higher  order  panel 
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method  from  the  constant  element  method  involved  the  derivation  of  new 
influence  coefficients  based  on  the  integration  of  a  variable  source 
density  over  a  curved  element.  Other  portions  of  the  method  were 
unchanged.  However,  the  development  of  the  higher  order  velocity 
equations  also  required  different  programming  logic. 

In  examining  the  potential  for  development  of  the  higher  order 
methods,  Hess  (1979)  noted  that  “a  consistent  approach  always  uses  a 
source  polynomial  one  degree  less  than  the  panel  polynomial."  Through  an 
independent  effort,  Brebbia  (1984)  presented  a  higher  order  approach 
using  the  direct  method  to  solve  for  a  surface  potential  polynomial 
stating  that  the  potential  function  must  be  of  a  degree  at  least  equal  to 
the  degree  of  the  polynomial  describing  the  element.  Knowing  that  the 
velocity  function  is  the  derivative  of  the  potential  function,  these  two 
observations  agree.  Asa  result  of  his  derivations,  Hess  showed  that  the 
solution  of  a  flat  element  with  a  constant  source  requires  one  integral,  a 
paraboloidal  panel  with  a  linearly  varying  source  density  requires  six 
integrals,  and  a  cubic  element  with  a  quadratic  source  density  requires 
twenty-three  integrals.  Development  of  higher  order  methods  has 
focused  on  the  paraboloidal  element  with  the  linearly  varying  surface,  as 
solutions  of  higher  order  than  that  offered  little  benefit  for  the  amount 
of  effort  required  to  produce  a  working  program  (Hess  1979).  Hess 
(1979)  and  Eriksson  (1983)  have  independently  developed  programs  for 
-  three  dimensional  higher  order  panel  methods.  The  higher  order  Hess 
program  evolved  from  the  constant  element  program  which  he  developed 
in  the  early  1960s,  while  Eriksson  developed  a  new  program  based  on  the 


work  of  Johnson  and  Rubber t  (1975).  Continued  work  in  the  near  future  is 
expected  to  deal  primarily  with  refinement  of  the  paraboloidal  element 
with  a  linearly  varying  source  (Eriksson  1983). 

In  order  to  alleviate  the  burden  of  preparing  the  input,  a  geometry 
package  for  input  data  generation  has  been  developed  which  is 
incorporated  into  the  Hess  higher  order  panel  program.  This  allows  the 
user  to  enter  relatively  few  points  to  describe  the  body.  The  geometry 
package  enhances  the  surface  representation  by  distributing  additional 
points  on  the  surface  based  on  one  of  many  algorithms  or  recurring 
geometries  (Halsey  1978). 

As  the  state  of  the  art  in  fluid  dynamics  has  progressed,  the  XYZ 
Potential  Flow  program  has  seen  increasingly  complex  applications 
requiring  a  great  deal  of  effort  in  preparing  the  input,  and  requiring  long 
computer  run  times.  Hess  (1979)  reported  the  use  of  the  Hess  constant 
element  program  for  a  configuration  utilizing  7000  effective  elements. 
Realizing  that  the  computation  time  increases  as  the  square  of  the  number 
of  elements,  it  is  easy  to  see  the  motivation  for  developing  the  higher 
order  panel  methods.  Though  modern  computers  offer  storage  capacities 
which  can  handle  most  applications  of  the  constant  element  method,  the 
higher  order  panel  methods  can  provide  equal  accuracy  for  much  less  user 
effort.  While  the  constant  element  method  is  still  a  versatile  tool, 

.  future  generations  of  the  surface  singularity  methods  will  be  able  to 
handle  the  more  complex  applications  being  demanded  in  fluid  dynamics. 


13.0  VELOCITY  CALCULATIONS  FOR  A  TRIAX1AL  ELLIPSOID 

As  the  only  true  body  for  which  an  analytical  solution  exists,  a 
triaxial  ellipsoid  was  selected  for  the  sample  calculations  in  order  to 
compare  calculated  results  with  the  analytical  solution.  Hess  has  made 
use  of  the  triaxial  ellipsoid  throughout  his  works  in  developing  both  the 
constant  element  method  and  the  higher  order  panel  method.  Therefore, 
the  XYZPF  program  will  be  compared  with  existing  results  of  the  Hess 
method  (Hess  1979). 

The  triaxial  ellipsoid  utilized  for  the  calculations  has  semiaxes 
dimensions  of  1.  2.  and  0.5  in  the  x.  y,  and  z  directions  respectively.  The 
surface  was  discretized  by  selecting  fixed  intervals  of  0.1  in  the  y 
direction,  and  fourteen  equal  divisions  of  the  90°  sector  in  the  x-z  plane. 
The  values  of  x  and  z  were  then  solved  in  terms  of  y  and  an  angle  0.  This 
method  yielded  280  elements  in  the  first  octant  fora  total  of  2240 
effective  elements  after  employing  symmetry.  AFORTRAN  program  was 
used  to  generate  the  corner  points  and  to  prepare  the  input  file  for  later 
use  by  the  XYZPF  program. 

Figures  (16)  and  (17)  show  excellent  correlation  with  the  analytical 
solution  and  little  difference  from  the  Hess  solution  using  4320  effective 
elements.  The  use  of  the  centroid  as  the  control  point  is  an 
approximation  used  to  simplify  the  multipole  expansion  of  the  potential 
about  the  origin  of  the  local  coordinate  system.  This  approximation  is 
valid  for  most  elements.  However,  for  elements  which  are  long  and  thin, 
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the  physical  difference  between  the  centroid  location  and  the  null  point 
location  is  significant,  and  use  of  the  centroid  can  produce  significant 
error  as  may  be  observed  in  figure  (16)  when  the  value  of  y  approaches 
2.0. 

Recent  calculations  on  the  same  body  (Hess  1979)  showed  that 
results  of  at  least  equal  accuracy  could  be  obtained  using  only  480 
effective  elements  using  the  higher  order  panel  method.  These  results  are 
a  significant  demonstration  of  the  value  of  the  higher  order  panel  method. 
Using  the  higher  order  panel  method  rather  than  the  constant  element 
method,  the  user  has  the  option  of  obtaining  equal  accuracy  with  cruder 
discretization  or  higher  accuracyfor  the  same  discretization  effort. 

While  the  results  of  the  triaxial  ellipsoid  show  relatively  little 
improvement  in  accuracy,  the  most  significant  advantages  of  the  higher 
order  panel  method  are  evident  for  a  body  with  concave  regions  (Hess 


Figure  16.  Comparison  of  analytic  and  calculated  velocity  distributions  on  an  ellipsoid  with 
axes  ratios  1  :2:0  5.  Velocities  in  the  xz-plane.  (from  Hess  and  Smith  1%2) 
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The  objectives  of  this  paper  were  (i)  to  describe  the  details  of  the 
approximation  of  an  arbitrary  three-dimensional  body  using  quadrilateral 
elements,  and  (2)  to  provide  a  detailed  derivation  of  the  exact  source 
panel  integrations.  Both  of  these  objectives  were  met. 

The  method  of  surface  discretization  and  source  panel  geometry  is 
easily  described  using  basic  principles  of  geometry  and  vector  algebra. 

By  using  quadrilateral  surface  elements,  many  surfaces  can  be  discretized 
in  a  very  straight  forward  logical  fashion.  The  user  can  frequently 
visualize  the  contour  lines  of  the  surface  which  may  be  used  to  form  the 
quadrilaterals,  with  some  help  from  an  intuitive  approach  to  the  fluid 
dynamics  problem.  The  method  of  forming  the  planar  quadrilateral 
element  in  the  XYZPF  Program  differs  slightly  from  the  method  presented 
by  Hess  and  Smith  (1962).  The  differences  lie  in  the  formation  of  the 
local  coordinate  system  and  the  use  of  the  centroid  rather  than  the  null 
point  as  the  control  point  for  applying  the  boundary  conditions.  The 
method  of  forming  the  local  coordinate  system  has  no  effect  on  the 
potential  flow  calculations  as  long  as  one  of  the  coordinate  vectors  is 
the  outer  normal  to  the  planar  element.  The  use  of  the  centroid  as  the 
control  point  is  an  approximation  used  to  simplify  the  multipole 
expansion  of  the  potential  about  the  origin  of  the  local  coordinate 
system.  This  approximation  is  valid  for  most  elements.  However,  for 
elements  which  are  long  and  thin,  the  physical  difference  between  the 
centroid  location  and  the  null  point  location  is  significant,  and  use  of  the 


centroid  can  produce  significant  error  as  may  be  observed  in  figure  (16) 
when  the  value  of  y  approaches  2.0. 

A  detailed  derivation  of  the  exact  source  panel  integration  has  not 
previously  appeared  in  literature,  though  the  results  are  summarized  by 
Hess  and  Smith  (1962).  The  derivations  presented  in  this  paper  verify  the 
equations  presented  by  Hess  and  Smith  (1962),  and  the  equations  used  in 
the  XYZPF  Program.  The  integral  expressions  for  the  velocity 
components  were  evaluated  exactly  with  no  assumptions  or 
approximations  used  in  the  course  of  the  integrations.  Since  the  method 
of  integration  reduces  the  surface  integral  to  a  line  integral  around  each 
of  the  sides  of  the  element,  the  integration  method  can  be  generalized  for 
a  planar  element  with  any  number  of  sides,  though  the  surface 
discretization  used  in  the  XYZPF  Program  uses  only  quadrilateral 
elements. 

The  calculation  of  potential  flow  about  arbitrary  three  dimensional 
bodies  is  an  engineering  tool  which  is  basic  to  design  involving  fluid 
dynamics.  The  XYZPF  Program  is  a  useful  tool  which  has  proven  its  value 
over  the  past  14  years.  However,  the  increasing  demands  placed  on  this 
method  are  exposing  the  errors  of  the  approximation  as  evident  in  the 
sample  calculations  presented  in  this  paper.  The  requirement  for 
increased  accuracy  has  motivated  the  development  of  the  higher  order 
_  panel  methods.  Some  of  the  limitations  imposed  on  the  XYZPF  Program 
were  due  to  computer  memory  and  speed  limitations.  Advancesin 
computer  performance  may  allow  future  investigators  to  eliminate  some 
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of  the  simplifying  approximations  used  in  the  XYZPF  Program,  allowing 
increased  accuracy  without  violating  computer  limitations.  Some 
modifications  might  include  the  use  of  the  null  point  as  the  control  point 
rather  than  the  panel  centroid  (as  is  used  In  the  Hess  program),  or 
extending  the  range  in  which  the  exact  velocity  calculations  are 
performed.  The  gains  in  accuracy  by  modifying  the  "constant  element 
method"  are  limited  by  the  basic  approximations  of  the  planar  element 
and  the  constant  source  density  for  each  element.  Significant  gains  are 
most  evident  in  the  higher-order  panel  methods.  This  author  concurs  with 
Eriksson  (1983)  in  expecting  advances  in  the  surface  singularity  methods 
■  o  focus  on  the  "development  and  refinement"  of  the  higher  order  panel 
methods. 
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APPEND 1 X  I  -  XVZPF  SECTION  PF 1 

PROGRAM  PFP  K I NPUT =128, OUTPUT  =128, TAPE5= I NPUT , TAPE6=0UTPUT , TAPEG3 , 
1  TRPE3=TRPE03, TRPE04, TRPE4=TAPE04, TAPE50= 128 ) 

XV2  POTENT  I RL  FLOW  PROGRAM  VERSION  4  SECTION  1 
READS  INPUT  AND  COMPUTES  QUADRILATERAL  PARAMETERS 

FOR  INFORMATION  CONTACT 

BILL  CHENG  OR  JANET  DEAN 
NUMERICAL  FLUID  DVNAMICS  BRANCH  CODE  1843 
NAUAL  SHIP  RESEARCH  AND  DEUELOPMENT  CENTER 
BETHESDR,  MARYLAND  20084 


DIMENSION  !NDEX<9,3),G(9,6),F<9),C2<9>,  IP(9>,XP<9>,VP<9)  ZP 

1  ,MSK< 100),  WS<24Q),PRQB< 15),DM(65Q> 

COMMON  X<800  ), Y<8GG  ),2<800  ), ID<41,?1 ),B<250),T<4600) 

EQUIVALENCE  (C2(1>,F<1>  ) 

EQU I  VALENCE < MS <  1 ),  KP<  1 ) ),  <MS<  10 1 ),  MSK<  1 )  >,  <WS<20 1 ),  NP  ),  (US 

2  NSP ),  <U$(203),NEP >,  (WSC204  >,NSE),  <US(205),MIX),  <US<206),IH' 

3  (MS  <20?  > ,  M I Z  >,  <US<209 ),  I  CUM ) ,  <148(2 10  >  I SM )  <US<2 11)  K)  (MS 

4  EPS >, (MS <208 >, IPS  ), <MS<2 12), I PF  >, <US(2 17), X I ), <WS<2 IS ) . V 
5(145(219), 2 1  ) 

EQU I  VALENCE  <  V 1 2, Y23 ) , < Y34 . V4 1 > 

INTEGER  P,P  1  ,P2,  P3,  F'4,PC  , P5,  P6 , P?,P8  ,P9 

WRITE  (6,5) 

5  FORMAT (49H 1XY2  POTENTIAL  FLOW  PROGRAM  SECTION  1  VERSION  4 
10  FORMAT  <11,  15A4  ) 

20  FORMAT  ( IX,  1514) 

50  FORMAT  <  1X,2I7,6E12.5) 

30  FORMAT  <1X,5F12.9> 

A.  READ  IN  CONTROL  PARAMETERS 
US <220  )=4 . 

K  1=0 
ID  1=0 
ID2=0 
103=0 
ID4=0 
I  D5=G 
1 06=0 
107=0 
MAXN=70 
MRXM=40 
MAXNQE=650 
MAXPC-800 
i CNTRL= 1 
E0F50=0 . 

READ  (5  ,  10)J,  (PROBd  ),  1  =  1,  15) 

IF  (EOF (5 )  ,EQ.  0. )  GO  TO  9 
WRITEC6,  8) 

8  FORMAT (39H0N0  TITLE  CARD  FOUND  -  PROGRAM  ABORTED  ) 

STOP 

9  CONTINUE 
J=0 

WRITE  (6,10)  J, (PROBd  ),  1  =  1, 15) 

SA=  .0 
SB=0 . 

21  FORMAT  (17H0N0.  OF  QUADS.  =,  14  /  17H  NO  OF  SEn 1 0NC;= 

2I4/31H  MAX.  NO  OF  ITERATIONS  X  FLOW  , I3,9H  Y  FLOW  "  ' 

313, 9H  2  FLOW, 13) 

READ  <5  , 1 1 )NQE,NSE,MIX,MIY,MI2, ISM, EPS, IUCT , IPS  IFF  ISP 


KP(IGG) 
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1  , 1 ED  I T 1 , I ED  I T3 , 1 ED  I T4 ,  I TftPE , XCENTER , VCENTER , 2CENTER 

11  F0RMATC6I4,F8.5,8I4, 1X,3F5.3> 

IF  (EOF (5)  .EQ.  0. )  GO  TO  19 
WRITEC6, 18) 

18  FORMAT C43H0N0  PARAMETER  CARD  FOUMD  -  PROGRAM  ABORTED  ) 

STOP 

19  CONTINUE 

IF  (IEDITJ.EQ. 1)  ICNTRL=0 
WRITE  <6, 2 1 )  NQE, NSE ,  M I X ,  M I V ,  M 1 2 
31  FORMAT  C'COHUERGENCE  CRITERIA",F8  5) 

41  FORMAT  (4H0  M,?X,2HX1 , 12X,2HX2, 12X.2HX3, 12X,2HX4, 12X,2HXP, 12X, 

1  2HXN, 12X, 1HA, 13X, 3HCZ4/4H  N,7X,2HY1, 12X,2HY2, 12X,2HV3, 12X,2HY4, 

2  12X,2HVP, 12X,2HYN, 12X,2HFL, 12X,3HCZ5/4H  P,?X,2H21, 12X,2HZ2, 12X, 

3  2H23. 12X,2HZ4,  12X,2HZF', 12X,2HZN, 12X,4HC21  , 10X,3HCZ6/> 

42  FORMAT  (4H1  M,7X,2HX1, 12X,2HX2, 12X,2HX3, 12X,2HX4,  12X.2HXP,  12X, 

1  2HX.N, 12X, 1HA, 13X,3HCZ4/4H  N,7X,2HV1, 12X,2HV2, 12X,2HY3, 12X,2HY4, 

2  12X , 2HVP . 12X . 2HYN , 12X . 2HFL .  12X, 3HC25/4H  P,?X, 2HZ 1 , 12X, 2HZ2, 12X, 

3  2HZ3, 12X,2HZ4, 12X,2HZP, 12X,2HZN, 12X,4HCZ1  , 10X,3HCZ6/> 

24  FORMAT  <  1H  , I1.19H  PLANES  OF  SYMMETRY) 

240  WRITE  (6,24)  ISM 
270  WRITE  (6.31)  EPS 
280  IF  (IPS.LE.O)  GO  TO  290 
285  WRITE  <6,36 >  IPS, IPF 

36  FORMAT  (45H0NEU  SOURCE  DENSITY  TO  BE  COMPUTED  FOR  QUADS., 1 4, 3H  -  , 
114) 

290  K=0 
WRITE(6J39)ISP 

39  FORMAT  <QHO  ISP  =  ,13  ) 

WRITE<6, 37)  I  ED  I T  1 , IEDIT3, IEDIT4, ITAPE 
37  FORMAT  (9H0IEDIT1  =,I3/9H  I  EDITS  =,I3/9H  IEDIT4  =,I3/9H  ITAPE  =  , 

1  13) 

WR I  TEC 6, 38 )  XCENTER . YCENTER . ZCENTER 
33  FORMAT (10HQXCENTER  =  .F5.2/10H  YCENTER  =,F5.2/10H  ZCENTER  =,F5.2) 
MM=0 
MN=0 
P=1 
0=1.0 

DO  291  1=1.41 
DO  291  J= 1 , 7 1 

291  I  DC  I . J)=0 
J=0 

C  B  READ  FIRST  PT . 

I ERR=0 

IF  Cl  TAPE . EQ . 1 )  GO  TO  292 
2000  READ  <5, 40)  XI ,YI . 2 1 ,NI ,MI , NS. ME, UN 

IF  CEOFC5).NE.O.  OF;'.  NS.LE.O)  GO  TO  2050 
WRITE (6,  45)  I CNTFiL , NS 
45  FORMAT C I  1 , 9H  SECT  I  ON  ,14; 

LINE=0 
GO  TO  293 

2050  IF< IERR.EO.O)  GO  TO  2200 
WR I TEC6, 2 100 ) 

2100  FORMAT C39H0N0  POINT  CAROS  FOUND  -  PROGRAM  RBORTED  ) 

STOP 

2200  I  ERR  =1 
I TRPE= 1 

WRITEC 6, 2300) 

2300  FORMAT C47H0ERR0R  IN  INPUT  -  POINT  CARDS  NOT  ON  INPUT  FILE  , 10X, 

1  53HPR0GRRM  WILL  CHANGE  ITAPE  TO  1  AND  TRY  TO  READ  TAPE50  ) 

IF  C  EOF <5 ) . EQ . 0  )  WRITEC6,2400)  XI,VI,ZI 

2400  FORMAT < 1 1HQEXTRR  FL0U,3F12.5,5X,2OHWILL  NOT  BE  COMPUTED  ) 
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292  READ (50, 40 >  XI ,VI ,ZI ,NI ,MI ,NS,NE,UN 
40  FORMAT  <3F12.9,4I4,F12.9> 

E0F50=E0F(50) 

IF  (E0F50.NE.0.  .OR.  NS.LE.O)  GO  TO  2450 
UR I TE<6, 45 )  ICNTRL,NS 
L I NE=0 
GO  TO  293 

2450  IF  (IERR.EQ.O)  GO  TO  2500 
UR  I TE(6 , 2 100  ) 

STOP 

2500  I ERR= 1 
I TAPE=0 
UR I TE (6, 2600) 

2600  FORMAT ("ERROR  IN  INPUT  -  POINT  CARDS  NOT  ON  TRPE50" ,  10X, 

V PROGRAM  UILL  CHANGE  ITAPE  TO  0  AND  TRV  TO  READ  INPUT  FILE") 

GO  TO  2000 

293  UNR=UN 
NSS=HS 
PC=1 

IF  (NE  . EQ . 0 )  GO  TO  2700 
Mf1 1  N=N  I 
MMAX=N I 
NM I N=M I 
NMRX=M 1 
GO  TO  300 
2700  MM  1 11= Ml 
MMRX=M I 
NM I N=N I 
NMRX=N I 
GO  TO  300 

295  IF  < I  TAPE . EO .  1 )  GO  TO  297 

READ  (5,40)  XI, VI, 21, Nl .Ml, NS, ME, W 
IF  (EOF(S).EQ.O. )  GO  TO  299 
NS=0 
Xl=0. 

V I  =0 . 

21=0. 

GO  TO  299 

297  READ(50,40>  Xi , VI ,21 ,NI .Ml , NS, ME, UN 
E0F5O=EOF (30 ) 

IF  (E0F50  NE.O)  NS=0 

299  PC=PC+1 

IF  (NS.HE.NSS)  GO  TO  330 

300  IF  (NE  EQ  0)  GO  TO  304 

301  I H=N I 
N I  =M  I 
M I  =  I U 

:  C.  STORE  PT.  IN  FT.  RRRRV 

304  IF  (MAXPC+  1-F'C )  295,305,310 

305  UR  I TE (6, 306)  NS,MI,NI 

306  FORMAT (60H  ERROR  IN  INPUT  -  THERE  ARE  TOO  MANV  DATA  POINTS  IN  SEC 
1 T I  ON  , 1 4, 30H  -  POINTS  BEGINNING  UITH  M  =, 14, 5H  N  =  ,14. 

2  17H  UILL  BE  IGNORED  ) 

L INE=L INE+ 1 
104=104+1 
GO  TO  295 
310  X(PC)=XI 
V(PC  )=V I 
Z(PC)=ZI 

IF  (MILE  MAXM  .AND.  NI.LE.MAXN)  GO  TO  315 
URITE(6,31t)  Ml ,NI 


LINE=LINE+  t 

311  FORMAT (38H  ERROR  IN  INPUT  -  I  HURL  ID  M,N  INDICES  , 10X, 

1  14HP0INT  WITH  N  =  , 14, 5H  N  =  ,14, 17H  WILL  BE  IGNORED  ) 

105=105+1 
PC=PC>  1 
GO  TO  295 
315  1 0 < M I ,  Nl >=PC 

MMAX=MAXO<MMAX,MI  > 

MMIN=MINO(MMIN,MI  ) 

NMfiX-MRXO<NMftX,NI  > 

NMIN=MINQ<NMIN,N! ) 

GO  TO  205 

330  IF  < I  ED  I T 1 . EG  1)  GO  TO  294 
IF  <LINE.LT.40>  GO  TO  333 
URITE<6,42> 

LINE=0 
GO  TO  294 

333  URITE<6,41> 

294  CONTINUE 

E.  00  LOOPS  TO  SWEEP  PT.  flRRftV 
N1=NNIN 

MM2=MMAX-MMIN 
NN2=NMRX-NM I N 

IF  (  MOL<MM2,2).EQ.O  . RND .  MOD<NN2.2).EQ.O)  GO  TO  332 

UR 1TE<6,331 )  NSS, MM  I N , MMflX , NM I N , NMftX 

LINE-LINE+1 

331  FORMAT ( 16H0ERR0R  -  SECTION  ,I5,45H  DOES  NOT  HRUE  QUADS  ARRANGED  IN 

1  BLOCKS  OF  4  ,9H  MM  I  N=  ,  I2,6H  MMAX=  , I2,6H  NM  I  N=  ,  12  , 

26H  NMftX®  ,12) 

I D6= I D6+ 1 

332  MM2-MM2/2 
NN2=NN2 }2 

DO  404  NN= 1 , NN2 
M  I=MM I H 

DO  402  MM= 1 , MM2 
NQ=1 

F.  HAUE  9  CORNER  PTS.  BEEN  GIUEN 

I T= I D<M 1 ,N 1 >* I D(M 1+ 1 . N 1  )* I D<H 1+2 . N 1 >* I D(M 1 . N 1+ 1 >* I D<H 1+ 1 , N 1+ 1  >* 

1  I D<M 1+ 1 , N 1+2  >* I D<M 1 , N1+2 >* ID<M 1+ 1 , N 1+2 >* 1 CKM 1+2,  N 1+2 > 

I FC  IT.EQ.O  )  GO  TO  402 
I  EF:R=0 
M2=M  1  + 1 

N2=N  1+ 1  =M  f ,  M2 
DO  400  N=N1,N2 
GO  TO  <334 ,335, 336,33?)  NC| 

334  P  1=1  CKM  ,N  ) 

P2*ID<M+1,N  > 

P3=ID(M+1,N+1> 

P4= I D(M  , N+ 1 ) 

P5=ID<M+2,N  > 

F'6=ID(M+2,N+1 ) 

P?=ID<M+1,N+2) 

PS=ID<M  . N+2 ) 

P9=P  1 

IF<<X<P1  ).NE  X<P2)  OR  V<P1  ).NE . VCP2).0R  2<P1 >.NE  2<P2>>  .AND 
1  <X<P1 ).NE.X(P4 ).0R. YCP1 ).NE . V<P4 ).0R  Z<P1  ).NE  Z<P4 )))  GO  TO  340 

F'9=IDCM+2,N+2> 

GO  TO  340 

335  P 1= I D<M  , N+ 1 ) 

P2= I D<M  ,  N  ) 

P3=ID<M+1 ,N  ) 


P4=  1 0(n+ 1 , N+ 1 ) 

P5=ID(M  H—  1  > 

P6=ID(M+1JN-1 > 

P7=ID(M+2,N  ) 

PS=ID(M+2,N+1 ) 

P9=P  1 

IFCCXCP 1 ) . NE . XCF'2 ) . OR  .  YCF'  1  > .NE .  YCF’2  > . OR .  ZCF'  1 ) . NE  .  ZCP2 ) )  .AND. 

1  (X(P1  ).NE.X(P4).0R.V(P1  ).HE.V(P4).0R.2(P1 >.NE.2(P4)>>  GO  TO  340 

P9= I DCM+2 , N—  1  > 

GO  TO  340 
P1=ID<f1+l,N  ) 

P2=ID01+1,N+1> 

P3= i D  <  M  ,  N+ 1 ) 

P4=IO(f1  ,N  ) 

P5=ID<M+1,N+2> 

P6=ID(M  ,N+2> 

P7=idcm-i,n+i) 

P8= I D<M— 1 , N  > 

F'9=P  f 

IF( <X<P1  >.HE .XCP2) .OR. V<P1 )  NE  V<P2).0R.2<P1 ).NE ,ZCF2>)  . RND . 

1  <X<F'  1  > .  NE . X<P4 ) . OR . V<P 1 ) . NE . V<P4 ) . OR . 2<P t  > . NE . 2<P4  > ) >  GO  TO  340 

P9= i D(H- 1 .N+2  > 

GO  TO  340 
’  P 1= I DCM+ 1 j  N+ 1 ) 

P2=!DCM  ’h+o 
P3=ID<M  ,  N  > 

P4=IDCM+1,N  ) 

P5=  I  DCH-  1 ,  H+ 1 ) 

pd=i:kh-i',h  > 

P7=IDCH  .  N- 1 ) 

F'8= !  D(M+  1 ,  N-  1  > 

P9=c'  1 

! FC CXCP f  ) . NE . XCF'2 ) . OR . VCP 1 > . NE . VCF2 > . OR . 2CP t  > . HE .  ZCF2 > >  . AND . 

1  CXCP  1  > . HE . XC.P4 > . OR . YCP 1 > . HE . V<P4 )  OR . 2(P  1 ) . HE  Z(P4 > )>  GO  TO  340 
P9= I D(H- 1 , H- 1 ) 

I  IPC1)=P1 
IP<2>=P2 
IPO  >=P3 
IF'(4)=P4 
IPO  )=P5 
IP(6)=F'6 
IPO  >=P7 
IP(8>=F'S 
I  PC 9  >=F'9 

02  COMPUTE  NORMAL  UECTOFi  <XN,YN,2N> 

X1=X(P3 >-X<P1> 

X2=X(P4 >-X(P2) 

Y1=VCP3>-V(P1> 

V2=VCF'4  >-VCF2  ) 

Z1=Z(F‘3>-Z<P1> 

Z2=ZCP4>-ZCP2> 

XH=V2fZ 1-V 1*22 
VN=X 1*Z2-X2*Z 1 
ZN=X2*V 1-X 1+V2 
R=S02CXH.YH,ZN> 

IF  CR.GT.  .00000000001)  GO  TO  345 
UR i TEC 6, 343) 

)  FORMAT C33H  ERROR  IN  INPUT  -  ZERO  AREA  QUAD  > 

L!NE=LINE+1 

101=101+1 
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xe=o. 

vc=o. 

2C=0 . 

FL=0 . 

CZC 1 )=0 . 

C2<4)=0. 

C2C5>=0. 

C2(6)=0. 

I ERR=  1 
GO  TO  351 
345  CONTINUE 
XN=XN/R 
YN=YN/F; 

ZN=ZN/R 
00=  5*R 

CONFUTE  CENTROID 
X1=XCP3  )-X(F'2) 

Y  1=YCP3 >-Y(P2> 

2 1=ZCP3  )-2(P2 ) 

X5=V  1*Z2-Y2*2 1 
Y5=Z 1*X2-22*X 1 
Z5=X1*Y2-X2*V1 
R1=S02<X5,Y5,Z5> 

R2=R-R 1 
IT=  1 

XC=  <  X  <  P2  >+X  <  P4  >+  <  ft  1  *X'  <  P3  >+fl2*X  <  P 1  >  >  /R  )  /3 . 
YC*WP2 >+VCP4 >+<fi  1*V(P3  •+F2+Y(D  1  >  >/F;  >/3 . 

ZC=  <  Z  <  P2  HZ  <  P4 )+  <  ft  1  +2  <  P3  >+fl2*  Z  <  P 1 ) )  /P, )  /3 . 
CONFUTE  SECOND  ftND  THIRD  UECTORS 
945  X4=VN+Z1-V1*ZN 
V4-ZN+X 1-2 1*XN 
Z4=XN*,v’1-X1*VN 
ft=  1 . /SQ2CX4 , Y4 , 24 ) 

X4=X4*fl 

Y4=V4*fl 

Z4=Z4*fi 

X3=ZN*V4-24*YN 
Y3=XN*Z4-X4*ZN 
23= V  N*  X4 -  V4  *  X'  N 

COMPUTE  POINTS  IN  QUftD  SYSTEM 
DO  94?  1=1,9 
L= I  PC  I  ) 

XFC I  )=X3*<X(L  )-XC:  >+V3* C VCL  >-YC  >+23*<Z<L  >-2C > 
VPC I  )=X4*<XCL )-XC )+V4*(VCL  )-VC >+Z4*(2<L >-ZC ) 
94?  ZPC I >=XN+CXCL >-XC )+YN+CVCL >-YC HZt1*<Z<L X-2C ) 
COMPUTE  MATRIX  COEF ,  TO  FIND  SURFACE  EQ. 

DO  943  1=2,9 
GCI,  0=1. 

GC 1 , 2  >=XF  C I ) 

G<  1 ,3>=VP( I  ) 

GC  1 , 4  )=XF'<  I  )**2 
GC  1 ,5)=VF'C  I  )**2 
GC I ,6>=YP< I >*XPC I ) 

949  F( I >=ZP< I ) 

DO  953  1=1,6 
GO,  I  )=GC9,  I  > 

GC5 , I  )=GC5, I )+GC6, 1 > 

953  GC6, I >=GC7, I >+GC8, I > 

FCO=F<9> 

F(5)=FC5)+F<6) 

F(6)=F(?)+F<8> 


'X'C'V-'X  ••• 
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SOLUE  MATRIX  EQ.  G*C2=F  FOR  C2 
CALL  MATINS<G,9,6,F,6, 1, DETERM, I  DM, INDEX) 

IF  (IDM.EQ.  I)  GO  TO  <955,960)  IT 
I ERR= 1 

URITE(6,954> 

954  FORMAT  <33H  ERROR  IN  INPUT  -  SINGULAR  MATRIX  ) 
L I N£=L I NE+ 1 

I D2= I D2+ 1 
GO  TO  060 

955  !T=2 

FIND  NEW  NORMAL  UECTOR 
XN=XN-C2<2 )*X3-C2<3 >*X4 
VN=VN-C2 <  2 )*Y3-CZ <  3 )*Y4 
2N=ZN-CZ  <  2 )*Z3-C2  <  3 )*24 
A=1./SQ2<XN,YN,2N) 

XN=XN*A 
YN=YN*A 
ZN=ZN*R 
GO  TO  945 

STORE  DATA 

960  EKJ+I  )=XP<  1 ) 

B<U+2 >=YF'( 1 ) 

B<J+3)=XP<2) 

B(J+4)=YP(2) 

B<J+5)=XP<3) 

B<  J+6  )=XP<4 ) 

B(J+7)=VP<4 ) 

B<J+8)=X3 
B<  J+9 )=Y3 
B(J+10)=23 
B< J+ 1 1 )=X4 
B< J+ 12 )=V4 
B<J+13)=24 
B<J+14)=CZ< 1 ) 
e<J+15)=CZ<4) 

B<J+ 16 )=C2<5 ) 

B<J+17)=CZ<6) 

IF  (K  LT.  7*MAXN0E )  GO  TO  965 
I D7= I D7+ 1 
K  1=K+K  1 
K=0 

965  CONTINUE 
T<K+1  )=XC 
T(K+2  )=VC 
T<K+3)=2C 
T(K+4  )=XN 
T<K+5  )=VN 
T  <K+6  )=ZN 
TCK+7 )=RQ 

COMPUTE  QUAORIJPOLE  MOMENTS 
XI  1=XP( 1 )+XP<2) 

X 1 2=XP< 1  )+XP<4 ) 

X 1 3=XP  <  3  >+XP ( 2 ) 

XI4=X'F'<3)+XP<4) 

XI5=XP<2)+XF'<4) 

VI  1=YP( 1 )+YP<2) 

Y!2=YP< 1 )+YP<4) 

YI3»VP<3)+YP<2) 

VI4=VP(3)+YP<4) 

V 1 5= VP  <  2  )+ VP  <  4  ) 

R1=A1/24. 


R2=fi2/24 

R3-AQ/12. 

AXX=(X I 1**2+X 1 2**2 >*R 1+<X 1 3**2+X 1 4**2  >*R2+X 1 5**2*R3 
fiXV=<X I 1*V I 1+X 1 2*V 1 2 )*R 1+<X  I 3*V 1 3+X 1 4*V 1 4  >*R2+X 1 5*V 1 5*R3 
AYY=< Y I 1**2+Y 1 2**2  >*R  1+<  Y I 3**2+ Y 1 4**2  >*R2+Y 1 5**2*R3 
COMPUTE  SOLID  ANGLE 
XX=XC-XCENTER 
YY=YC-YCENTER 
22=2C-ZCENTER 
X  1«XX*X3+YY*Y3+ZZ*Z3 
V  1=XX*X4+YY*Y4+22*Z4 
2 1=XX*XN+VV*YN+22*2N 
RD= 1 . /S02<X 1 ,  Y 1 ,  Z 1 ) 

RCU=RD**3 

RSU=RCU**2*RD 

Sfl=SA+Z 1 *  <  RO*RCU- ( ( AXX* ( Y 1 **2+2 1**2-4 . *X 1**2  ) 

1  +RYY*<X 1**2+Z 1**2-4 . *Y 1**2 >  >*  1 . 5- 15 .  *X  1*Y 1*RXY )*RSU ) 

e(.J+ 18  >=AXX 
EKJ+19)=AXV 
B< J+20 )=RVY 
ERROR  TESTS 

01=SQ2<«P<3)-XP<  1  >), (YP(3)-VP<  1  >),0.  ) 

D2=SQ2<<XP<4 )-XP<2 > ), < YP<4 )-VP<2 )), 0 . > 

FL= . 5+RMRX 1  (D 1 , 02  > 

CZ23=ABS<  CZ< 2)  )  +  RBS<  CZ<3>  ) 

IF<  RBS < CZ  <  2 ) )+RBS <  CZ  <  3 ) )  .GT.  FL*.001  >  GO  TO  970 


IF  (  ABS(  CZ<1» 


FL*. 3  )  GO  TO  97? 


975  FORMAT <29H  QUESTIONABLE  POINT  -POOR  FIT  .6E14.3) 
I  EFlR=  1 
L I NE=L I NE+ 1 

977  IF  (YP< 1 )*XP<2 >-VP<2)*XP< 1 )  .GE.  0.  . RND . 

1  VP<2)*XP<3>-YP<3>*XF'<2>  -GE.  0.  .AND. 

2  YP<3)*XP<4>-YP<4)*XP<3)  .GE.  0.  .AND. 


VP  <  4  >*XP  <  1  >- VP  <  1  >+XP  (4  >  .  GE .  0 . 


GO  TO  984 


980  WRITE <6,  1000)  <  XP<  I  >,  VP<  t  >,  1  =  1,4) 

1000  F0RHRT<41H  ERROR  IN  INPUT  -  CROSSED  OR  CONCAUE  QUAD 
1  4<2F  10 . 5, 3X>  > 

I ERR=  1 
L I NE=L I NE+ 1 
1 03= 1 03+ 1 

984  CRCF-SQ2C <XP<2 )-XP< 1 ) ) , <VP<2 >-YP< 1 > ), 0  >+XP<3 )-XP<2 >+  SQ2< <XP< 1 > 
1  XP <4»,< VP <  1  )-VP < 4 > > , 0 .  >+SQ2 < < XP < 4 )-XP (3>>,< VP < 4 >- VP < 3 ) ) , 0 

IF  <  36 . *RQ  GT.  CRCF**2  )  GO  TO  986 
LINE=LINE+1 
WRITE*.  6, 95 1 ) 

961  FORMAT  <24H  WARNING  LONG  THIN  QUAD.) 

986  IF  (  Z1  .GE.  0.  )  GO  TO  351 

350  WRITE  <6,85) 

85  FORMAT  <35H  QUESTIONABLE  POINT  -  INWARD  NORMAL) 

I  ERF;=  1 
L I NE=L I NE+ 1 

J.  EDITE  QUAD  INFORM AT 'ON 

351  IF  < I  ED  I T 1 . EQ . 2 . AND .  IERR.EQ  0)  GO  TO  354 
IF  <  I  ED ! T 1  EQ.  1  )  GO  TO  354 

GO  TO  (356,357,358,359)  NO 

356  WRITE<6,51)  M,X<P1).X<P2).X(P3),X<P4).XC.XN.ftQ  ,CZ<4)  , 


GO  TO  360 
35?  URITE<6, 5 1 ) 


M,  X(P 1 ),X(P2),X(P3 ) ,X(P4 ) , X C . XN , AQ  , CZ<4 ) 

N,  Y<P  1  >, Y<P2), Y<P3>, Y<P4 ), VC, YN.FL  ,CZ<5) 
P, ZCP 1 ), Z(P2 ),Z(P3),Z(P4 ), ZC, ZN, CZ< 1 ), CZ<6 ) 

M, X(P2 ), X<P3 ), X<P4 ), X<P 1 ) . XC . XN , RQ  ,CZ<4) 


1 

2 


N,V<P2),V<P3),V<P4),V<P1),VC,VN,FL  ,02(5) 

P,Z<P2),Z<P3  ),Z<P4  ),Z<P1  ),2C, ZN,CZ<  1  ),CZ<6 ) 


GO  TO  360 

358  WRITE(6,51 )  M,X<P4),X<P1 ),X<P2),X<P3),XC,XN,RQ  ,CZ<4) 

1  N ,  Y(F'4  >,  Y<P  1  >,  V<P2  > .  Y(P3  > , YC . YN .  FL  ,CZ<5) 

2  P,Z<P4),Z(P1)/Z(P2),Z(P3)/ZC;ZN;cZ<1>;CZ(6) 
GO  TO  360 

359  WRITE<6,51)  f1,X<P3>,X<P4>,X<F1  >,X<P2),XC,XN,flQ  ,CZ<4) 

1  N,V<P3),Y<P4>,Y<P1>.Y<P2>,YC,YN.FL  ,CZ<5) 

2  P,Z<P3),Z<P4),Z<P1>,Z<P2),ZC,ZN,CZ< 1>,CZ<6> 

360  CONTINUE 

51  FORMAT  <1H  , I3,8E14.5/1X, I3,8E14.5/1X, l3,8E14.5/> 
LINE=LINE+4 

IF  (LINE.LT. 50)  00  TO  354 
352  WRITE  <6, 42) 

LINE=0 

354  CONTINUE 
JaJ+20 
l=P 

DM< I  )=UNR 
P=F'+ 1 
NC!=NQ+  1 
IERR=Q 
349  K=K+7 

K.  WRITE  OUT  BLOCK  OF  B  RRRfiV  IF  FULL 
IF  <J.LT .240)  GO  TO  400 

355  WRITE  (04)  0, <B< I >, 1=1,240) 

Q=P 

J=0 

L.  END  OF  DO  LOOP  OUER  PT.  RRRfiV 
400  CONTINUE 

402  M  1=M 1+2 

404  N  1=N 1+2 

M  SET  FOR  NEXT  SECTION 
NSS=NS 

DO  405  M=MMIN,MMfiX 
DO  405  N=HM!N,Nf1fiX 

405  ID<ri,N)=0 
F'C=  1 

IF  <ME  EO  0)  GO  TO  410 
MMfiX=N I 
MM  I  fi-N  I 
NM I N=M I 
NMft <=M I 
GO  TO  420 
410  MMfiX=MI 
MM  I N=M I 
NM I N=N I 
NMflX=N I 
420  NE=ME 
UNFl=UN 

IF  (NS.LE  0)  GO  TO  500 
WRITE <6, 45)  ICNTRL.NS 
LIHE=Ci 
GO  TO  300 

500  WRITE  T4)  Q,  <B<  I  ),  1  =  1,240) 

550  NP=<K+K 1 )/? 

I SMP  =  ISM  +  1 

GO  TO  <595,590,580,570), I SMP 
570  $fi=Sfi+Sfi 
580  Sfi=Sfi+Sfi 


590  SA=SA+SA 

01  WRITE  PARAMETERS  AND  T  ARRAV  ON  TAPE 
595  J  =  1 

IF  (ITAPE.EQ.  1  AND.  E0F50  .NE.O  >  GO  TO  601 
59?  WS(J)  =  XI 
USCJ+20)  =  VI 
WS( J+40  >  -  21 
J  =  J+1 

IF(XI**2  +  Vl**2  +  21**2)  599,599,598 

598  WRITE(6,600)  XI, VI, 21 

600  FORMAT ( 1 1H0EXTRA  FLOW, 10X,3F  12.5) 

601  READ(5,40)  XI, VI, 21 

IF  <E0F<5>  ,EQ.  0)  GO  TO  597 
Xl=0. 

Vl=0. 

21=0. 

GO  TO  59? 

599  IF  < I SP . LT . 0 >  GO  TO  605 
WRITE  (03)  (PROBd  ),  1  =  1,  15) 

WRITE  (03)  (WSd ), 1=1,220). IE0IT3, IEQIT4 
WRITE  (03)  (Td  ),  1  =  1, K> 

WRITE  (03)  (DMd  ),  I  =  1 , NP ) 

605  CONTINUE 

N1  CHECK  SOLID  ANGLE 
6 1u  WRITt  (6, SO)  SA 
80  FORMAT ( 1 4H0S0L I D  ANGLE  =  ,F8.3) 

620  REWIND  04 
REWIND  03 

622  I DS= I D 1 + 1 D2+ 1 03+ 1 D4+ 1 D5+ 1 D6+ 1 D? 

IF  (IDS.E0.0  AND.  NF.EQ.NOE)  GO  TO  638 
WRITE (6, 625) 

IF  ( I D 1  . GT . 0 )  WRITE(6,628)  101 

IF  ( ID2  . GT . 0 )  URITE(6,629)  ID2 
IF  ( ID3  . GT . 0 )  WRITE(6,630)  ID3 
IF  ( ID4  GT.  0)  URITE(6,631 )  ID4 

IF  (IDS  .GT.  0)  WF;ITE(6,632)  IDS 

IF  ( 1 06  .GT.  0)  WF;ITE(6,633)  ID6 

IF  (ID?  .GT.  0)  WRITE (6, 634)  NF' , MAXNOE 

IF  (NP.NE.NQE)  WRITE (6, 637)  NF' , NQE 
STOP 

625  FORMAT (3SH0FATAL  ERROR  IN  DATA  -  PROGRAM  ABORTED) 

628  FORMAT ( 1H0, 1 5, 3 1H  QUAOR I  LATERALS  W I TH  ZERO  AREA  > 

629  FORMRT( 1H0, I5,44H  QUADRILATERALS  GENERATE  A  SINGULAR  MATRIX  ) 

630  FORMAT (1 HO, IS, 26H  CROSSED  QUADR I  LATERALS  ) 

631  FORMAT ( 1 HO, 15, 32H  SECTIONS  HAVE  TOO  MANY  POINTS  ) 

632  FORMAT ( 1H0, 15, 34H  PO I  NTS  HRUE  I NUAL ID  M.N  I  HD  I CES  ) 

633  FORMAT ( 1H0, IS, 52H  SECTIONS  DO  NOT  HRUE  OUADS  ARRANGED  IN  GROUPS  0 

IF  4  ) 

634  FORMAT ( 1H0, I5.4SH  OURDR I  LATERALS  GIUEN,  EXCEEDING  THE  LIMIT  OF  , 

1  14) 

637  FORMAT  ( 1  HO ,  1 5 . 28H  OUADF;  I  LATERALS  G I UEN ,  NOT  ,14) 

638  IF  (I SP  LE . 0 >  GO  TO  640 
WRITE (6, 639)  ISP 

639  FORMAT (7H0  ISP=  , 14,  19H  -  PROGRAM  ABORTED  ) 

STOP 

640  CONTINUE 

02  READ  PEPS2  AND  TRANSFER  TO  IT 
STOP  1 
END 


FUNCTION  SQ2(X, V,Z) 


ooo  oooooonoooooooooooo 


C  COMPUTE  SQUflR  ROOT  OF  R**2 

R*  ABS(X)+ABS(V>  +ABS(Z>  +.0000000000001 
700  RS=X**2+Y**2  +2**2 
R=R+RS/R 
R= . 23+R+RS/R 
R=R+RS/R 
SQ2=.25*R+RS/R 
RETURN 
END 

SUBROUT INE  MRT I  NS ( A , NR , N 1 , B , NC , M 1 , DETERM , ID, INDEX) 

MATRIX  INUERSION  UITH  ACCOMPANVING  SOLUTION  OF  SIMUL.  EQ. 
PIUOT  METHOD 

FORTRAN  IU  SINGLE  PRECISION  UITH  ADJUSTABLE  DIMENSION 
NOVEMBER  1971  S  GOOD  NAUAL  SHIP  R  &  D  CENTER 
UHERE  CALLING  PROGRAM  MUST  INCLUDE 

DIMENSION  A(NR,NR>,  B(NR,NC),  INDEX(NR,3) 

UHERE  NR,NC  ARE  DIMENSIONS  OF  A, B, INDEX 
N 1  IS  THE  ORDER  OF  A 

Ml  IS  THE  NUMBER  OF  COLUMN  VECTORS  IN  B  (MR','  BE  0) 
DETERM  UILL  CONTAIN  DETERMINANT  ON  EXIT 
ID  UILL  BE  SET  BV  ROUTINE  TO  2  IF  MATRIX  A  IS 
SINGULAR,  1  IF  INUERSION  UfiS  SUCCESSFUL 
MATRIX  R  (INPUT  MATRIX)  UILL  BE  REPLACED  BV  A  INU 
MATRIX  B:  THE  COLUMN  UECTORS  UILL  BE  REPLACED 
BV  CORRESPONDING  SOLUTION  UECTORS 
INDEX'  UORKING  STORAGE  ARRAV 

IF  IT  IS  DESIRED  TO  SCALE.  THE  DETERMINANT  CARD  29  MAV  BE 
DELETED  AND  DETERM  PRESET  BEFORE  ENTERING  THE  ROUTINE 

DIMENSION  A(NR,NR>,  B(NR,NC>,  INDEX(NR.S) 

EQUIVALENCE  (IROU,JROU),  ( ICOLUM, JCOLUM),  (RMAX,  T,  SUAP) 

INITIALIZATION 

N=N  1 

n=m 

DETERM*  1.0 
DO  20  J=  1 ,  H 

20  I NDEX( J , 3 )=0 
DO  550  I  *  1 ,  N 
C 

C  SEARCH  FOR  PIUOT  ELEMENT 
C: 

RMAX  =0.0 
DO  105  J= 1 , N 

IF(.INDEX(J,3)-1 )  60,  105,  60 
60  DO  100  K=1,N 

IF(iN0EX(K,3>-1>  80,  100,  715 
80  IF  (  AMR)'  -ABS  (A(J,K)))  85,  100,  100 
85  IROU  =  J 
ICOLUM  =  K 
RMAX  =  ABS  (R(J,K)> 

100  CONTINUE 
105  CONTINUE 

INDEX (ICOLUM,  3)  =  INDEX( ICOLUM, 3)  +  1 
INDEX  1,1)-  IROU 
INDEX (I, 2)  =  ICOLUM 
C 

C  INTERCHANGE  ROUS  TO  PUT  PIUOT  ELEMENT  ON  DIAGONAL 


v r»v^_ 


i 

i 


c 


E 


I' 

■ 


C 

c 

c 


c 

c 

c 


i 


c 


I 


c 


IF  < I ROW- I COLUM )  140,  310,  140 
140  OETERM=  -DETERI1 
DO  200  L=1,N 
SUftP=  ft( IROU,L) 

RC  IROU,L  )=RC  I  COLUM,  L  > 

200  R< I COLUM, L )=SWflP 
IFCM)  310,  310,  210 
210  DO  250  L= 1 , M 
SWflP-B(IROW,L) 

BC IROU,L)=BC I COLUM, L) 

250  BC I COLUM, L )=SUAP 

010 IDE  P 1 00 T  ROW  BV  PIOOT  ELEMENT 

310  PIOOT  =  RC  ICOLUfl,  I  COLUM) 

DETERM=OETERM*P  I OOT 
330  RCI COLUM,  ICOLUfl)  *  1.0 
DO  350  1=1, N 

350  RC I COLUM, L >=RC I COLUM, L )/P I OOT 
IF  CM)  380,  380,  360 
360  DO  370  L=l' M 

370  B  C I COLUM , L )=B< I COLUM , L ) /P I OOT 

REDUCE  NON-PI OOT  ROUS 

380  DO  550  L 1= 1 , N 

I FCL 1-1 COLUM)  400,  550,  400 
400  T=RCL 1 , I COLUM > 

RCL1, I COLUM >=0.0 
DO  450  L= 1 ,N 

450  P C L 1 , L >=A<L  1 , L  )-RC I COLUM, L  >*T 
IF  CM)  550,  550,  460 
460  DO  500  L=  1 , M 

500  B  C L 1 , L  >=B  C L 1 , L  >-B ( I COLUM , L  >*T 
550  CONTINUE 

INTERCHANGE  COLUMNS 

DO  710  1=1, N 
L=N+ 1 - 1 

IF  < I NDEX <L, 1 >-l NDEX C L, 2 > >  630 ,  710,  630 
630  UR0U=  I NDEX CL,  1) 

JCOLUM= I NDEX CL , 2 ) 

DO  705  K= 1 , N 

SURF1  =  RCK,  JFlOU  > 

ft  C  K , JROW >=R C  K , JCOLUM > 

RCK, JCOLUM >=SURF 
705  CONTINUE 
710  CONTINUE 
DO  730  K= 1 , N 

I FC I NDEXCK ,  3 )- 1  >  715  ,  720  ,  715 
720  CONTINUE 
730  CONTINUE 
ID- 1 

810  RETURN 
715  ID=2 


GO  TO  810 


« 


I 


PROGRAM  PF2(G(JTPUT= 128, TAPE6=0UTPUT, TRPE03 , T RPE02, TRPEOS, 

1  TAPE09 ,  TRPE0 1 ,  TAPE  1 1 ,  TAPE04 ,  TRPE4=TRF'E04 , 

3  TAPE 1 -TAPED 1 , 

2  TAPE3-TAPE03,  TAPE2-TAPE02, TAPE8-TAPE08,  TAPE9-TAPE09  > 

C 

C  XV2  POTEMTIRL  FLOW  PROGRAM  VERSION  4  SECTION  2 
C  COMPUTES  MATRIX  COEFFICIENTS 
C 

COMMON  B(241),T(46G0>,  UK  1000), V2(  1000), U3(  1000), CK  900), C2(  9 
200),C3(  900),  WS<300  ),PROB( 15 ), 

3UX(8),VY(S).VZ<8> 

EQIJ I  VALENCE ( Y3 ,  V2  >  ,  < US < 20 1  > ,  NP  > ,  ( US ( 2 1  CO,  SVM ) ,  ( US < 2 1 1 ) .  KM  > 

1  ,  (WSC212).  IPF)  ,  ( WS  (  208  >,  I  PS  >  ,  <  KM ,  MK  > 

INTEGER  SVM 
BLK=  1 . 0 
1014-0 

READC03)  (PRGB(I ), 1=1, 15) 

UR  I TE(6, 5 ) 

5  FORMAT (49H0XV2  POTENTIAL  FLOW  PROGRAM  SECTION  2,  VERSION  4  ) 

90  FORMAT ( 1H0, 15A4) 

WRITE  <6  ,90)  (PRQB<  I >. 1  =  1, 15) 

C  A.  READ  PARAMETERS,  T  ARRAY,  FIRST  BLOCK  OF  B  ARRAY 

READ(03)  (USC I ), 1  =  1,220) 

RE AO (03)  (T< I >,  1-1 , MK) 

READ (04)  (8(1 >, I = 1 ‘ 24 1> 

C  B.  START  LOOP  OVER  QUADRALATERALS 

K=  1 
■J=1 
P-1 
•JC=  1 
JV=2 

NPN-5+<<NP+4>/5) 

KMM-7+NFN+ 1 

290  I F(B(  1  )-P  )595 , 295,595 

98  FORMAT ( " PO I  NTS  OUT  OF  ORDER  B(  1  )=",  1F4.0, "  P= " ,  1F4 . 0 ) 

295  J*2 
290  XI-B(J) 

V  1=E'(  J+  I ) 

X2=E( J+2 ) 

V2=p'  ,.!♦?,  > 

X3=B( J+4 > 

C  V3=V2 

X4-BCJ+5 > 

V4=B(  J+6  '• 

XN=T(K*3> 

YN=T (K+4  ■ 

2N=T(K+5 ) 

XC=T’>  ■ 

VC-T  .;k+  1  > 

2C=T(K+2 ) 

A=T (K+C ) 

XX-B'..  J+7  1 
VX«B< J+8 ) 

ZX=B( J+9 ) 

X'V=B(.J+ 10 ) 

VV=6v J+  1 1  > 

2Y=B(J+12) 

C  Cl  COMPUTE  LENGTH  OF  SIDES  OF  QUAD 

0 1 2=?02F ( X 1 ,  X2,  V  1 ,  Y2 . 0  ,0  ) 

023-SQ2F ( X2 . X3 . V2 . V3 , . 0 . . 0 ) 


C34=SQ2F < X3 , X4 ,  V3 ,  V4 ,  .0,  .0) 

D41«SQ2F<X4,X1,V4,V1,  .0, .0) 

C2  COMPUTE  SLOPE  OF  SIDES 
I F<X2-X3 >305,300, 305 

300  0123=1. 

GO  TO  310 

305  CM23=  <  V2-V3  > / < X2-X3  > 

Cl  23=0. 

310  IFCX3-X4 >315,31 1,315 

311  C 1 34= 1 . 

GO  TO  320 

315  Ct134=<V4-V3 )/(X4-X3 > 

C 1 34=0 . 

320  I FCX4-X  1)325,321,325 

321  CI41=1. 

GO  TO  330 

325  CM4 1=< V 1— V4 )/(X  1-X4 ) 

C 14 1=0 . 

330  IFCX1-X2  >335 ,331, 335 

331  Cl  12=1. 

GO  TO  340 

335  CM  1 2=  C  V2- V 1 > / <  X2-X 1  > 

Cl  12=0. 

C3  COMPUTE  QURDRAPOLE  MOMENTS 

340  C I XX=B( J+ 17  > 

CIXV=B(U+18> 

CIVV=B(J+ 19 > 

CV 12=0.0 
CXI 2=0.0 
CV23=0 ■ 0 
CX23=G . 0 
CV34=0  0 
CX34=0 . 0 
CV4 1=0 . o 
CX4 1=0.0 

C4  COMPUTE  SIM  AMD  COS  OF  SLOPE  ANGLE  FOR  EACH  SIDE 
IF  CD  12 >  9341,9342,9341 

9341  CV12=(V2-V 1 )/D 12 
CX 12=CX 1— X2  >/D 12 

9342  I F ( D23 >9343 , 9344 , 9343 

9343  CV23=(V3-V2 >/D23 
CX23=  < X2-X3  >/D23 

9344  I F  < 034  >9345 , 9345 ,9345 

9345  CV34= < Y4-V3  > /D34 
CX34=<X3-X4 >/D34 

9345  I F ( 04 1 >9347 , 934? / 9347 

9347  CV4 1  =  < V 1 - V4  > /D4 1 " 

CX4 1=(X4-X 1 >/D4 1 

C5  COMPUTE  MAX  LENGTH  OF  QUAD 

9343  ST =ABS  <  X 1 -X3  > 

ST 2-SQ2F (. X2 ,  X4 ,  V2 ,  V4 ,  .0,  C  > 

ST=RMAX  1  <  ST , ST2 , 0 1 2 , 023 . D34 ,041) 

0  START  LOOP  0"ER  NULL  PTS 

342  K0= 1 
L=  1 

343  I =K0 

I F  ( I  F’S  >  9350 , 9350 , 9350 

9350  IF(.L-IPS>  9305,9352,9352 


9352  IFCL-IPF)  9350,9350,9355 
9355  C 1 : JC )= . 0 
C2  ( UC  >= .  0 
C3( JC >=  0 


GO  TO  541 
IS=  1 

XCQ=T ( I > 

VCQ=T  <  I  + 1 ) 

2C0=T< 1+2) 

XNQ=T ( I  +3 ) 

VNQ=T  <  I  +4 ) 

ZNQ=T <  I  +5  ) 

E.  COMPUTE  DISTANCE  BETWEEN  QUAD  AND  NULL  PT. 

.  DETERM  IN  METHOD 
RPQ=SQ2F <XC,XCQ,VC, VCQ ,20,200) 

IF<RPQ-ST*4 )35G, 350,460 
X=«CQ-XC )*XX+(YCQ-VC )+YX+ <200-20 )*ZX 
V= ( XCQ-XC )*XV+ ( YCQ- VC )*V V+ (ZCQ-ZC >*ZV 
2=(XC0-XC >*XN+<VCQ-VC >*VN+< ZCQ-ZC >*2N 
I F <RP0-ST*2 . 0 )355  355,400 

F  COMPUTE  UELOC I  TV  COEF .  BV  EXACT  METHOD 


R  1=S02F (X, XI, V, VI, 2.0.  ) 
R2=SQ2F  <  X , X2 , V , V2 ' 2 , 0 .  ) 
R3=SQ2F <  X , X3 , V , V3, 2 , 0 . ) 
R4=S02F <  X , X4 , V , V4 ,2,0.  ) 
IF  ((R1+R2)  .LE.  D 12 >  i 


GO  TO 
GO  TO 


GO  TO 
GO  TO 


IF  <<R2+R3)  LE.  D23)  GO  TO  1000 

IF  <<R3+R4)  .LE.  D34)  GO  TO  1000 

IF  <CR4+R1>  LE.  04 1)  GO  TO  1000 

OLA  1  =RLOO <  <  R 1 +R2-D 12 )  / <  R 1 +R2+D 1 2  >  > 

CLA2=AL0&< (R2+R2-D23 )/(R2+R3+023 ) ) 

CLA3=RL0G( CR3+R4-Q34 )/<R3+R4+D34 ) ) 

CLR4=RL0G<  <R4+Fi  1-04  1  >/<F;4+R  1+04  1  >  > 

TUX=CV  12*CLfi  1+CV23*CLA2+CV34*CLA3+CV4 1+CLA4 
TUV=CX 12*CLA1+CX23*CLA2+CX34*CLA3+CX4 1+CLA4 
TUZ-0 . 

I F  <  ABS  <  2 /ST )- . 0 1 0  >  375 , 36 1,361 
230=2**2 

E  1=ZS0+(X-X 1 )**2 
E2=2S0+<X-X2 )**2 
E3=2S0+ (X-X3 )**2 
E4=ZSQ+(X-X4  )**2 
H  1  =  <V-V  1  >*<X-X 1 ) 

H2=  <  V-V2 )*<X-X2) 

H3=  < V-V3 >*<X-X3 ) 

H4*<V-V4 >*<X~X4 > 

IFCCi  12)363,353,364 
US1  =  <C'M12*E  1-H1  )/<Z*F;1 ) 

U  32=  <  CM  1 2+E2-H2 )  /  <  Z+Fi2 ) 

AT  1=ftTRN(WS 1 > 

AT2=ATAN<WS2> 

TUZ=AT 1-RT2 
IFCCI 23  >365 , 366 , 36? 

AT3=ATAN< CCM23+E2-H2 )/<Z*R2 ) ) 

AT4=RTRN< <CM23*E3-H3 )/<2*R3 > ) 

TvZ=TUZ+AT3-AT4 
IF  CO  1 34)363, 363,360 
BT5=ATAN<  CCM34+E3-H3 >/<2*R3 ) ) 

AT6=ATRN< (CN34+E4-H4 )/<Z+R4 ) > 

TUZ=TUZ+RT5-RT6 
IF CO  1 4 1)370 ,370,375 
AT7=ATAM<  COM4 1+E4-H4  )/<Z*R4  ) ) 

AT8=ATANC  COM4 1*E 1-H 1 )/(Z*R 1 ) ) 

TU2=TU2+RT7-AT8 
GO  TO  450 

G.  COMPUTE  UELOC I  TV  COEF.  BV  QUAORfiPOLE  METHOD 


RPQ3=RPQ**3 

RPQ7=<RPQ3**2>*RPQ 

WS1=  X/RPQ3 

X$Q=X**2 

YS0=Y*+2 

ZSQ=Z**2 

PS=YSQ+ZSC'-4 .  *XSQ 
QS=XSQ+ZSQ-4 . *YSQ 
US2=X* t 9  *PS+30 . *XSQ ) /RPQ7 
US3=3 .  *Y*PS/RPQ7 
WS4=3 . *X*0S/RPQ7 

TVX=R*US  1-C  I XY+US3-C I  XX+US2-C  I  YV*US4 
WS  1=Y/RP03 

US2=Y*(9.  *05+30 .  *VSQ  >  /RPQ7 

TW=ft*WS  1-C  I XX+US3-C  I  XY*US4-C  I YV+US2 

TUZ=Z*  C  R  /RFQ3-3 .  *  (  C I XX+FS-5  .  *C  I XY+X+Y+C I YY+QS  )  /'RPC7  > 

UX< I S )=TUX+XX+TUY+XY+TV2  *XN 

UV< IS  )=TUX*YX+TVY+YY+TV2*VH 

UZ( ! S )=TVX*ZX+TVY*ZY+TUZ*ZN 

GO  TO  470 

H  COMPUTE  VELOCITY  COEF.  BY  MONOPOLE  METHOD 
RRPQ3=fi / ( RPQ**3 ) 

VX <  I S  >=  < XCO-XC  X+firiF  Q3 

VY< IS)=  (VCO-VC )+RRPQ3 
UZC I S  i-  C2CQ-ZC )+ARFQ3 

I.  REFLECT  NULL  PT.  IN  PLANE  OF  SYMETRY 
GO  T0<480, 485, 490,495., 500,505,5 10,5 15 ),  IS 

DO  LOOF'8  SET  UF  TO  FORCE  USE  OF  INDEX  REGISTERS 
01=00 
J2=JC 

i  inv=ux(  i  > 

UD Y= Y Y (  1 ) 

UDZ=UZ i  1 ) 

VKJ1  >=VX(  1,> 

U2< J 1  )=UX( 1 > 

U3CJ1  )=VX< 1 ) 

V  1(01+1  )=UV( 1 ) 

U2CJ1+1 )=VY< 1 ) 

V3U1+1  >=VY<.  1  > 

U  1(0 1+2 )=UZ( 1 ) 

U2(  J 1+2  )=|JZ<  1 ) 

V3 < 01+2 )=UZ< 1 > 

IF(SVM)  530 , 530,48 1 
IS=2 

XZ  SYMETRY 

VCQ=-VCQ 
GO  TO  345 

I F ( SYM- 1 >5 1 7 , 5 1 7 , 485 
XV  SYMETRY 

I  S=3 

ZCQ=-ZCQ 
GO  TO  345 
1 8=4 

YCQ=-VCQ 
GO  TO  345 

I FCSVM-2 >5 16,5 15 . 435 
VZ  SYMETRY 

I  S=5 

XCQ=-XCQ 
GO  TO  345 
IS=6 

VCQ=-VCQ 


GO  TO  345 
505  IS=7 
2CQ=-2CQ 
GO  TO  345 
510  !S=8 
VCQ=-YCQ 
GO  TO  345 

J.  ROD  CONTRIBUTIONS  OF  RLL  REFLECTIONS 

515  U  K  J 1  >=U  K  J 1  HUX  <S HUX ( 7 HUX ( 6 )+UX ( 5 ) 

02(01 >-02 ( J 1 >-UX<8 HUX  <7 >+0X<6 HUX  <5 > 

03 v J 1 >=03 (J 1 HUX(8HUX<7 HUX(6 >+0X(5 > 

0 K J 1+ 1 >=U 1 ( J 1+ 1 HUVC8  HUV<7  )+0V<6  HUVC5 > 

02  <  J 1  + 1  >=02  (01+1  HUY  ( 8  >+UV  (  7  HUY  <  6  HUY  (  5  > 

03  (J  1+ 1  >=03  ( J 1+ 1  HOY  <8  HOY  <7  HUY  <6  HUY  <5  > 

U 1  <  J 1  +2  >=U  K  J 1 +2  >-02  <  8  HUZ  <  7  HU2  ( 6  >+U2  <  5  > 

02  ( J 1  +2  >=02  ( J 1  +2  HUZ  <  8  HUZ  <  7  HUZ  ( 6  HUZ  <  5  > 

U3<  J 1+2  >=U3(.J  1+2  HUZ  <8  HUZ  <7  HUZ  <6  HU2C5) 

516  UK  J 1  >=U  K  J 1  HUX  <  4  HUX  (  3  > 

02  <  J 1  >=02  ( J 1  >+0X  <  4  HUX  (  3  > 

03(J1 >=U3<J1 )-UX(4  )-UX<3 > 

U  1(01+1  >=0 1  (  J 1+ 1  HUY  ( 4  >-UY  ( 3  > 

02  ( J 1  + 1  >=02  <  J 1  + 1  HUY  ( 4  HUY  <  3  > 

U3  <  J 1  +  1  >-03  <  J 1+  1  HUY  <  4  HUY  <  3  > 

0  K  J 1 +2  >=0 1  ( J 1+2  HUZ  <  4  HUZ  <  3  > 

02 (J  1+2  >=U2C J 1+2  HUZC4 HUZ<3 > 

03  ( .J  1  +2  >=03  <  J 1  +2  HUZ  < 4  HUZ  <  3  > 

517  0 1 ( J 1 >=U 1 (J 1 >+0X(2  > 

U2(J1  )=U2(J1 >-0X<2> 

03  <  J 1  >=03(01  HUX  <2> 

0  K  J 1+ 1  >=0  K  J 1+ 1  HOY  <2  > 

02<  J1+ 1  >=U2<  J1+ 1  HUY  <2  > 

03  ( J 1  +  1  >=03  ( J 1  + 1  HOY  ( 2  > 

Kite  1+2  >=U  1 C J 1+2  >+UZ(2  > 

02  ( J 1  +2  >=U2  ( J 1 +2  HUZ  <  2  > 

03  ( J 1 +2  >=03  <  J 1 +2  HUZ  <  2  > 

530  C 1 <J2 >=XNQ*U 1< J 1 HVNO+U 1 C J 1+ 1 >+ZNQ+0 1 C J 1+2 > 

C2 < J2 >=XNQ*U2 < J 1 >+YNQ*U2 <  J 1 + 1 >+ZN0*U2 < J 1 +2  > 

C3CJ2  >=XNC!+03(J  1  >+VNQ*U3<J1+ 1  >+ZNQ+03(J1+2  > 

540  00=00+3 
54  1  JC=J C+ 1 

0  UR ! TE  COEF I C ! ENTS 

K.  WRITE  COEF.  ON  TRPE  OR  DRUM  IF  STOP, ROE  AREA  IS  FULL 
545  !  F  >,  00- 1 00 1  >570 ,555, 555 

555  00=2 

OKI  >=E:LK. 

02 ( 1  >=BLK 
03 < 1 >=BLK 

I F  < BLK-636 . 0  >  560,553, 566 
553  WriiTE  <01>  ELK, 01, 02, 03 
GO  TO  568 
562  REWIND  01 

566  URI TE<  1 1 >  ELK. 01, 02,03 
568  BLK=BLK+ 1 . 

570  IF <0090 1)580,571,571 

571  I DW= I DU+ 1 

WRITE  < 02  >  I  DU , C  1 
WRITE  (08 >  IDW.C2 
WRITE  (00 >  I  DU, C3 
576  JC=  1 
580  K0=KQ+7 
L=L+ 1 

L  END  OF  LOOP  OUER  NULL  PTS 


IFCKQ-k'M  >243  .,343,531 
58)  CKJC>=0 
C-2<  JC  >=0 
C3<JC>=0 

IFCKQ-KMM >541 ,585,585 

585  P=F'+ 1 
K=k'+7 
J=J+20 

1FCK-KM  >536,586,600 

N.  END  OF  LOOP  ODER  QUADS, 

H.  READ  NEXT  BLOCK  OF  B  RRRRV  IF  NEEDED 

586  I F<J-241  >206,500,590 
590  READ' 04  >(B< I > , I  =  1 , 24 1 > 

J=2 

IF(B< 1 >-P >395, 296 .595 
595  WRITE  (6,98)  B<  1  >,P 
STOP 

600  IFCBLK-636 .0 >  610,620,630 

O.  WRITE  RENA  INI  NO  COEF,  ON  TAPE  OR  DRUM 
610  WRITEC01 >BLK,U1,U2,U3 

REWIND  01 
00  TO  640 
620  REWIND  01 

630  URITE< 1 1 >  BLK,U 1,02,03 
REWIND  11 

640  WRITE  (02 >  IDW.C1 

j  «r-, «  rr  i  <-•'  •  r--“- 

wn  i  IL  \  L'O  .•*  1  L.'M  J  Ua. 

URITE  <09 >  IDW.C3 
REWIND  02 
REM ! ND  03 
REWIND  04 
REWIND  OS 
REWIND  09 

P.  TRANSFER  TO  PFPS 3 
GO  TO  5000 

1000  WRITE (6, 2000)  L,P 

2000  FORMAT (3H  L=  . I5.20X.3H  P=  ,F5.1> 

5000  CONT I  HUE 
STOP  2 
END 

FUNCT I  ON  SQ2F (XI, X2 , V 1 , V2 , 2 1 , 22  > 

X=X 1-X2 


R$=Z**2+V**2+X**2 

R=RBSO<  >+ABS(V >+ABS(Z >+  1 , 0E-20 

R=R+F;S/R 

R= .  25*F;+F!S/R 

R=  R+RS/R 

SQ2F=  , 25*R+RS/P 


o  o  o  o 


PROGRAM  PFP3(0UTPUT= 128, 7RPEQ2, TRPE08, TRPEOQ, 

1  TAPE  1 2 , TRPE03 , TRPE6=0UTPUT , TRPE2=TAPE02 , 

2  TAPES=TAPE08 , TRPE9=TRPEQ9 , TRPE3=TRPE03 ) 

XV2  POTENTIAL  FLOW  PROGRAM  VERSION  1  SECTION  3 
SOLOES  MATRIX  EQUATION  FOR  SOURCE  DENS] TV 

COMMON  SN<654>,VIP<650>,S(5,650>,PRQB( 15>,WS(220),DM<65D>, 

1  B  < 220 ) , COEF  <  900  > , XN < 650 ) , VN < 650 ) , ZN  < 650) 

EQU i UALENCE  <US(213),EPS>  , <KK,B<201)> 

EOU I UALENCE  <MIX,WS<205)>, (Ml Y,WS<206>>, <mZ,US<20?)> 

EQU  I  UALENCE  (WSC201  >.NP  ),  043(208),  IPS),  <WS<212),  IPF),  CI4SC21 1  >,KM>, 
KKM.MK) 

5  FORMAT (49H1XVZ  POTENT  I RL  FLOW  PROGRAM  SECTION  3,  VERSION  4  ) 

WRITE  (6,5) 

READ  (03XPR0B<I  ),  1  =  1, 15) 

WR I TE  (  6, 100 1 )(PROB( I ), I = 1 , 15 ) 

1001  FORMAT <1 HO, 15A4) 

READ  (03 )  <M$< I  ), 1  =  1 , 220 ), I EO I T3, I ED  I T4 

READ  (03XSKIP, SKIP, SKIP, XN< I  ),VN( I  ),ZN( I  ),SK IP, I  =  1  .NP ) 

D=- .5/3. 14159265 

READ  <03)(DM( I ), 1=1, NP) 

K 1=  1 
K2=NP 

240  FORMAT  < 18H0CHANGES  IN  PROB  -, 15A4) 

IF( IPS)  1220, 1220, 1231 
123 1  READ  ( 12 X  BOO, K=  1,15) 

WRITE  <  6,240 ><B(K ), K= 1 , 15 > 

READ  < 12 )<  B(K),K=  1,220) 

READ  ( 1 2  )SK I F 

READ  <12 '.'SKIP 

K  1=1  PS 
K2= I PF 

C  A.  SET  CONDITIONS  FOR  FLOW  OF  -1  IN  X  DIRECTION 

1220  FX=- 1 
FV=0 
FZ=0 
NF=- 1 

C  B.  COMPUTE  INITIAL  APPROXIMATION  TO  THE  SOURCE 

1240  DO  1250  K= 1 , NP 

V I P(K  )=XN(K)*FX+YN(K )+FV+2N(K )*FZ-DM(K ) 

S(5,K)=-UIP(K>*. 11936 

C  C.  SET  PARTIAL  SUM  VECTOR  TO  ZERO 

1250  SN(K)=0. 

SN(NP+ 1 >=0 . 

SN(NP+2  )=0 . 

SN(NP+3)=0 . 

SN(NP+4  )=0 . 

WRITE <6, 997)  FX,FV,FZ 
URITE  (6,998) 

99S  FORMAT (27H0 1 TERftT ! ON  SUM  OF  CHANGES  , 9X ,  1HA , 10X , 2HB 1 , 10X , 2HB 2 ) 

IT=  1 
I  C=5 

IF  (IPS)  1260,1260,1255 

1255  READ(  12)  (S(5,K),  K=1,KK) 

DO  1256  K- 1 , KK 

DO  1256  1=1,4 

1256  S(l ,K)=8(5,K) 

C  D.  START  ITERATION 

1260  BANB=0 


IF  CNF)  1261, 1262,  1263 

1261  READ  (02)IDW,C0EF 
GO  TO  1264 

1262  READ  (08)IDU,C0EF 
GO  TO  1264 

1263  READ  <09)IDU,C0EF 

1264  J=0 

D.  READ  FIRST  BLOCK  OF  COEF 

E.  START  LOOP  OUER  QUADS. 

DO  1290  K= 1 , NP 

F.  PICK  UP  SOURCE  DENS  I  TV 
SP=S<IC,IO 

G.  START  LOOP  OUER  NULL  PTS. 

DO  1290  KP=1,NP,5 

I FC 0-900)80,65,65 
65  IF  CNF >67,68,69 
6?  READ  (02)  I  DU. COEF 
GO  TO  70 

6S  READ  (08) I  DU, COEF 
GO  TO  70 

69  READ  ( 09)  I  DU.  COC 

70  J=0 

H.  COMPUTE  PARTIAL.  SUM  FOR  NEXT  5  I 
80  SNCKF  )=SN(KP )+COEF( J+ 1 >*SP 

SNCKP+ 1 )=SNCKP+ 1 )+COEFC J+2 )*SP 
SNCKP+2 )=SNCKP+2 )+C0EF( J+3  )*SP 
SNCKP+3 )=SN<KP+3 )+COEF( J+4 >*SP 
SNCKP+4 )=SN(KP+4 HCOEFC J+5 )*SP 
J=J+5 

J.  END  OF  LOOP  OUER  NULL  PTS. 

K.  END  OF  LOOP  OUER  QUADS. 

1290  CONTINUE 

L.  COMPUTE  NEW  SOURCE 
IF  CNF)  91,92,93 

91  REWIND  02 
GO  TO  94 

92  REWIND  08 
GO  TO  94 

93  REMIND  09 

94  PRSS= 1 . 0 
SUM=0 . 

DO  100  K=K  1 , K2 
SNCK>=<  SNCK)+UIPCK)  )*D 
TEST=ABSC£NCK  )-SC I C, K ) ) 

SUt1=S'JM+TEST 

IF  (TEST  ,GT.  EPS)  PASS=-1.0 
100  CONTINUE 

IF  (PASS  .EQ.  1.0)  GO  TO  180 
IF  (IT.GE.MIX)  GO  TO  ISO 
IF  (IEDIT3  .EO.  0)  WRITEC6.99)  IT, SUM 
!T= I T+ 1 
I C= I C- 1 

IF  CIC  .EQ.  0)  GO  TO  120 
00  110  K=K 1 , K2 
S<IC,K)«SN<K) 

110  SNCK)=0. 

GO  TO  1260 
120  A=0 


DO  140  K=  K1,K2 
0S9=2*S< 1 , K )-SN(K  >-S(2, K ) 

IF(DS9  .GT.  0. )  GO  TO  122 
fi=A+S(2JIO-S(1,K> 

0R=DA-DS9 
GO  TO  125 

ft  =A  +S<1,K)-S(2,K> 

DA=DA+DS9 

0S1-S(4,K)-S(3,K) 

DS2=S<3,  K >-S<2, K ) 

DS3=DS 1-0S2 
DSS=S(2,K)-S(1,K) 

DS5=DS2-DSS 

DS6=DS 1-QSS 

0$4=DS2-SC 1,K)+SN(IO 

D$?=DS3*0S4-DS5*DS6 

DS8=DS6*D85-DS4*DS3 

IF<DS7  .GT.  0. )  GO  TO  128 

B  1=E  1-DS 1*DS4+DS2*DS6 

D1=D1-DS7 

GO  TO  130 

B1=B1+0S1*DS4-DS2*DS6 

D1=D1+DS7 

IF  (OSS  . GT.  0. )  GO  TO  132 
B2=B2-DS 1*DS5+D$2*DS3 
D2=D2-DSS 
GO  TO  140 

B2-B2+DS 1 *0S5-DS2*DS3 

D2=D2+DS8 

CONTINUE 

R=R/DR 

E'1=B1/D1 

B2=B2/D2 

IF( IT  ,EQ.  6)  GO  TO  155 

Hfi=H-fiS 

AA=ABS(AA> 

IF  CRli  .GT.  . 02 >  GO  TO  148 
DO  145  K=K 1 , K2 

S(5.K)=fl*(SN(K)-S<  I.KJHSC  1,10 
SN(K)=0. 

WRITE <6, 6000) 

FORMAT <  29X , 1 7Hfl  EXTRflPOLAT I  ON  ) 

GO  TO  160 
BB  1=B  1-BS 1 
BB1=ftBS(BB1 > 

BB 1 =50  *BB  1 

BB2=B2-BS2 

EB2=fleS(BE2) 

BB2=50 . *BB2 
BBB=RBS(B 1 >  +  RBS(B2) 

IF  (  (BB  1  .GT.  BBB )  OF; .  (BB2  .GT.  BBB )  >  GO  TO  1 
DO  150  K=K  1 , K2 

S ( 5 . K >=S ( 2 , K  )*B  1*<$(  t,K >-S < 2 , K > >+B2* <  SH < k >-S <2 ,  K > > 
SH(io=0. 

WRITE C6, 7000) 

FORMAT (29X, 17HB  EXTRAPOLATION  ) 

00  TO  160 
DO  158  K=K  1 , K2 
8(5,  K )=SN(K  > 

SN(K)=0. 

1 0=5 


URITEC6,  161)  ft,  Bt,  B2 
161  FORMAT <29X, 3E 12.3) 
flS=fl 
BS1=B1 
BS2=B2 
GO  TO  1260 

180  URITE<6,99>  IT, SUM 
DO  182  K=K  1 , K2 
182  S<  1,K)=SN(K> 

MR  I  TE<03 )  <S<  1,K>,  K-1,NP> 

99  FORMAT <4X, I3,E18.5) 

907  FORMAT  <  13H0  X  UELOC ITV=,F4  1, 15H  V  UEL0CiTV=,F4  1, 
1 15H  2  UELOC !TV=,F4  1) 

I FXFZ  >140D  1390,  1400 

PI  IF  THIS  UAS  NOT  LAST  FLOW,  SET  FOFl  NEXT  FLOW 
1390  FZ=FY 
FY=FX 
FX=0 
M I X=M I V 

Mm  V  —  j  i  i  *_ 

NF=NF+ 1 
GO  TO  1240 
1400  REMIND  03 

P2  REftD  IN  PFPS4  AND  TRANSFER  TO  IT 


o  o  o  o  o 


APPEND  I X  10  -  IVZPF  SECT  I  OH  PF4 

PROGRAM  PFP4 < OUTPUT , TRPE6=0UTPUT , TRPE03 , TRPE0 1 , TAPE  1 1 , 
1  TAPE3=TAPE03 , TRPE 1 =TAPEO 1 ) 

XY2  POTENTIAL  FLOW  PROGRAM  UERSION  1  SECTION  4 
COMPUTES  UELOCITIES  AND  PRESSURE  COEFFICIENTS  FOR 
POINTS  ON  THE  BODV 


COMMON  UX 1 <650 ),UV1 (650 >,UZ 1(650 >,  UX2<650>,W2(650>,UZ2<650> 

1  , UX3 < 650 ) , UY3 ( 650 ) , UZ3 ( 650 ) ,  SK650),  82(650),  S3C650) 

2  ,  X<650),  V<650),  Z<650),  T4(650),  T5(650),  T6(650) 

3  ,  DM ( 650  > 

D I  MENS  I  ON  PROS ( 1 5 ) , US ( 220 ) , CU 1 ( 1 000  > , CU2 ( 1 000  > , CU3 ( 1 000 ) 

EQU I UALENCE  (WS(201  ),NP),  (US(21 1  >,KM>,  OJS(217),UXI  ),  (WS(218),UVI  ) 
1(148(219), UZI  ), (I4SC208),  IPS),  (US(212),  IPF) 

EQU I UALENCE  (M IX,WS(205 ) ), (Ml V, US (206 ) ), (M I Z, US(20? ) ) 

5  FORMAT (49H0XVZ  POTENTIAL  FLOW  PROGRAM  SECTION  4,  UERSION  4  ) 

WRITE  (6,5) 

RERD  (OSXPROBd  ),  1  =  1,  15) 

100  FORMftT ( 1H  , 15R4 ) 

WRITE  (  6, 100)(PR0B(I >, 1=1,15) 

A  READ  PARAMETERS  AND  SOURCE 
READ  (03)  (US( I  ), 1=1,220), IEDIT3, IEDIT4 

READ  <03)<X<  I  ),V<  I  ),2(  I  ),T4(  i  ),T5(  I  ),T6(  I  ),SKIP,  1  =  1,  NP) 

READ  (03)(DM(I >, l=1,NP) 

D=-  5/3. 14159265 

READ  (03XSKI  ),  1  =  1, NP) 

READ  (03X82(1  >,  1  =  1,  NP) 

READ  (03)  (S3< 1  ), l  =  1,NP) 

J=1 
K 1=  1 
K2=NP 

IF< I  PS >108, 108, 102 
102  K 1= ; F  S 
K2= I PF 
108  BBFl=  1 . 

E.  READ  FIRST  BLOCK  OF  COEF. 

ITAPE=01 

RERD  (01)  BB,CU1,CU2,CU3 
IF  (BB-BBFl )  300,  120,300 
120  DO  125I-K1.K2 


UX  1  ( 1  )=- 1 . 0 

-SKI  >*T4(  1  ) 

/D 

UV  1  ( 

1  )= 

-SKI  )*T5(I  ) 

/D 

UZK 

1  )= 

-SKI  >*T6<  1 ) 

/D 

UX2( 

1  )= 

-S2( 1  >*T4( 1  ) 

/D 

UV2( 

1  )=- 1 . 0 

-S2( 1 >*T5( 1  ) 

/D 

UZ2< 

1  )= 

-S2( 1 )*T6( 1  ) 

/D 

UX3( 

1  )= 

-S3( 1 )*T4( 1  ) 

/D 

UV3< 

1  )= 

-S3( 1  >*T5< 1  ) 

/D 

125  UZ3( 

1  >*-1.0 

— S3( 1 )*T6( 1  ) 

/O 

C. 

SET  UP 

LOOP  OUER  QUADS. 

JC=2 

D. 

PICK  UP 

SOURCE 

130  S1J=SKJ) 

S2J=S2(J) 

S3J=S3(J) 

E.  SET  UP  LOOP  OUER  NULL  PTS. 

00  180  JP=K 1 , K2 

F  COMPUTE  PARTIAL  SUM  FOR  3  COMPONENTS  OF  3  UELOCITIES 
UX 1 ( UP )=UX K  JP  >+S  1 J*CU 1 ( JC ) 

UV 1  ( JP  )=UY  1  ( JP  )+S  1  J*CU  KJC+1 ) 
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U2 1 (JP  )=U2  K  JP  )+S  1  J*CU  KJC+2 > 

UX2( JP  )=UX2( JP >+S2J*CU2( JC ) 

UV2( JP )=UY2< JP >+S2J*CU2( JC+ 1 ) 

UZ2 ( JP )=UZ2 (  JP )+S2 J*  CU2<JC+2) 

UX3  < JP )= UX3 ( JP >+S3 J*CU3  <  JC ) 

UY3CJP )=UY3( JP )+S3J*CU3( JC+ 1  > 

U23<  JP  >=UZ3( JP  >+S3J*CU3( JC+2 ) 

JC=JC+3 

C  G.  READ  MORE  COEF.  IF  NEEDED. 

IF  (JC- 1000 >180, 135, 135 
135  JC=2 

IF(B6R-635.0)  140,150,155 
140  READ  (01)  BB,CU1,CU2,CU3 
GO  TO  160 
150  REWIND  01 

155  READ  (11)  BE, CU 1 , CU2, CU3 
160  BBR=BBR+ 1 . 

IF  (BBR-BB )  300, 180,300 
C  H.  END  OF  LOOP  ODER  NULL  PTS. 

180  CONTINUE 
J= J+ 1 

C  I .  END  OF  LOOP  OUER  QUADS . 

IF(J-NP)13Q, 130,200 
200  IF(BBR-635.0)  231,231,232 

231  REWIND  01 
GO  TO  233 

232  REWIND  11 

C  K.  EDIT  THE  UELOCITIES  ETC.  AND  WRITE  THEM  ON  TAPE 

233  WRITE  (03 )(UX 1(1 ), UY  1(  I ) , UZ 1  ( I )  ,|=1.NP) 

WR I TE  (03  KUX2(  I  ), UV2(  I  ),  U22(  I  )  ,  I  =  1 ,  NP ) 

WRITE  (G3)(UX3(I ),UY3(I ),U23(I )  ,I=1,NP) 

235  FORMAT ( 1H1, 15A4,8H  PAGE  =,115) 

REWIND  03 
IP=K1+49 
IS=K  1 
I PRGE= 1 

IF  ( IEDIT4  . EO.  1)  GO  TO  203 
IF  (MIX  . LE.O  )  GO  TO  265 
242  FORMAT (8H0  X  FLOW) 

240  FORMAT (4H  PT.,  10X, 1HX,9X, 1HV,9X,  1HZ.  DX , 2HUX, 8.X, 2HUV, 8X . 2HUZ . 10X, 
1  5HAB3.U,  ?X, 2HCP , 6X , 6HS0URCE , 4X , 8HU  NORMAL) 

245  FORMAT  ( IX, ! 3, 4X , 3F 10  5, 4X, 3F 10 . 5, IF  13 . 5, 2F 1 1  5.E12.2) 

250  I F ( I P-K2 >255 , 255 , 260 

C  J.  COMPUTE  PRESSURE  AND  ABS.  UALUE  OF  UELOCITV 

255  WRITE  (6, 235 >(PROB( I ), 1  =  1, 15), IPAGE 

WR i TE  (6,242) 

WRITE  (6,240) 

DO  257  l=IS, IP 

USQ=UX  1  ( I  )**2+UV 1 ( I )++2+UZ 1 ( I )**2 

UM=  ( ABS  ( UX  1(1)  )+ABS  ( UV  1(1)  >+ABS  ( UZ  1(D)  >* .  79 

UM=UM+USQ /UM 

UM=  25+UM+USQ/UH 

UM  = .  5*  ( UM+USO  /U11 ) 

CP  =  1 . -USQ 

UNR  =UX  1  ( I  )*T4( I )  +UY 1  ( I )*T5( I  )  +UZ 1  ( I  >*T6( I  ) 

257  WRITE  (6,245)  I ,X( I ),Y( I ),Z< I ),UX1( I ),UV1( I ),UZ1( I ),UM 
1  ,  CP  , SKI), UNR 

1 3= I S+50 
I P= I  P+50 
I PRGE= I PAGE+ 1 
IFCK2-IS)  265,260,250 
260  IP=K2 
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GO  TO  255 
265  IP-K1+49 
I  S=K  1 

IF  (Ml V  .LE  O  )  GO  TO  280 
26?  t  F  < I P-K2  >2?5 ,  273 , 2?0 
270  IP=K2 

275  WRITE  <6,235><PR0B< I >, 1  =  1, 15),  I PAGE 

WRITE  <6,27?) 

27?  FORMAT (SHO  V  FLOW) 

WRITE  <6,240) 

DO  278  l=IS, IP 

USQ=UX2< I )*+2+UV2< I >**2  +UZ2< I >**2 

UM=  <FIBS<UX2< I )>+ABS<UV2< I >>+RBS<UZ2<l >>>*  79 

UM=UM+U$0  /UM 

UM= .  25+UM+USQ/UM 

UM  =.5*<UM+US0/UM) 

UNR  =UX2<  I  )*T4<  I )  +UV2<  I  >*T5<  I  )  +U22<  I  >*T6<  I  ) 

CP  *  1 . -USO 

278  WRITE  <6,245)  I ,X< I ),V< I  ),2< I ),UX2< I ),UV2< I >,UZ2< I ),UM 
1  ,  CF  ,S2< I >,UHR 

I S= I S+50 
I P= I  P+50 
I PflGE= I PRGE+  1 
I F  < I S-K2 >267,267,280 
280  ! F-K 1+49 
I  S=K  1 

IF  <MI2  LE.O  >  GO  TO  293 
282  I F< I P-K2 >290,290,285 
285  I P=K2 

290  WRITE  <6,235 XPFlOEK  I  ),  1  =  1,  15),  IPRGE 

WRITE  <6,292) 

292  FORMAT  OHO  2  FLOW) 

WRITE  <6,240) 

00  291  (=18, (P 

US0=UX3< I >**2+UV3< I )**2+UZ3( I )**2 

UM=<RES<UX3< I  ))+RBS<UV3<l ))+RBS<UZ3<l >>>*.79 

IJH=UM+USC!/UM 

UM= . 25*UM+USQA!M 

UM  = . 5*<UM+USQ/UM ) 

CP  =  1 . -USQ 

UMR  =UX3 < I  )*T4< I )  +UV3 < I  >*T5  < ! )  +UZ3  <  t  )+T6< I ) 

291  WRITE  <6,245)  I , X < I ) , V < I > , 2 < I ) , UX3 < I >, UV3 < I ) , UZ3 < I ), UM 

1  ,CP  , S3( I  ) , UMR 

18=18+50 

IF- I P+50 

I  F'AGE=  I PRUE+ 1 

IF< I S-K2  >282 , 282, 293 

L.  CHECK  FOR  fl  FOURTH  FLOW 

293  J  =  1 

294  UXI  =  WSCJ) 

UV!  =  WS< J+20 ) 

UZI  =  WS< J+4G ) 

295  IF  <UX I **2+W I **2+UZ 1 **2 >  400,400,301 
301  I S=K 1 

IP=K 1+49 

320  IF  C I P-K2 >330, 325 , 325 

N .  EDIT  THE  UELOCITV  RMD  PRESSURE  FOR  FOURTH  FLOW 
325  I P=K2 

330  WRITE  <6,235)<PR0B<I ), 1=1, 15), IPRC€ 

WRITE  <6,315  >UX i , UV I , UZ I 

WRITE  <6,340) 

M.  COMPUTE  UELOCITV  AMD  PRESSURE  FOR  FOURTH  FLOW 


DO  333  l=IS, IP 

0X4  —  (OX I *OX 1 ( I HOY I *0X2 ( I >+02 1 *0X3  ( I  >  > 

0V4  =-  (OX  I  *0V  1  ( I  HUY  I  *0V2 ( I  >+02 1 *0Y3 ( I >> 

024  =-  (OX  I *02 1 ( I  HOY I *0Z2( I HOZ I *0Z3( I >  > 

030=0X4  **2+0Y4  **2+024  **2 

0M=(flBS(UX4  )+RBS(0V4  >+RBS(0Z4  >>*.79 

Of1=OM+OSO/OM 

0M= . 25+0M+0SQ/0M 

OM  = . 5* ( UM+USQ /UM ) 

CP  =  1 ,-(OSQ>/(  OY I **2+02 1 **2+0X I **2  > 
i  WRITE  (6,345>  I  > X ( I  >,V(  I  >,2(  I  >,0X4  ,074  ,024  ,011 

1  ,CP 
1 3= 1 3+50 
I P= I  P+50 
I PflGE= I PRGE+ 1 
I F ( K2- 1 S  >350 , 325 , 320 
)  J  =  J+1 
GO  TO  294 

i  FORMRTC 19H0  ONSET  FLOW,  UXI=,F6.3,2X,4H0VI=,F6.3,2X,4HUZI=,F6.3 
I  F0Rf1ftT<4H  PT., tOX, 1HX,9X, IHY,9X, 1H2, t3X,2H0X,8X,2H0V,8X,2H02,  WX 
1  5HRBS.0,  7X , 2HCP > 

i  FORMAT  (IX, I3,4X,3F  10. 5,4X,3F 10.5,  IF  13. 5, 1FJ1.5) 

I  CONTINUE 
GO  TO  5000 
1  CONTINUE 

WRITE  <6,310>  ITRPE,BBR,BB 

)  FORMAT  (6H1TRPE  , 12,  17H  OUT  OF  POSSITION/1 14.6F6.  1  > 

I  CONTINUE 
STOP  4 
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PROGRAM  PFP5C I  NPUT=  1 28 ,  0UTPUT=  1 28 ,  TRPE03 ,  TRPE04 , 

1TflF'E5=  I HPUT ,  TRPE6=0UTPUT ,  TAPE3=TRPE03  , TAPE4=TAPE04 ) 

XV2  POTENT  I RL  FLOW  PROGRAM  UERSION  4  SECTION  5 
COMPUTES  VELOCITIES  AND  PRESSURE  COEFFICIENTS  FOR 
OFF  BODY  POINTS 

COMMON  B  < 24 1  > ,  XP  <  500  ),YP< 500 ) . ZP  <  500  > , UX 1 <  500 ) , US  < 220 ) 

1 ,  VV 1  <  500  >,U2K 500  > ,  V.X2  ( 500 > ,  UV2  < 500  ) ,  U22  ( 500 ) ,  UX3  ( 500  > ,  UV3  (500) 

2,  UZ3C500  ),  UX<8) ,  UV<8  > ,  UZ<S),S  1  <650 > .  S2C650 ),  S3  <650 ),  V 1  <3 ),  V2<3  >,  U3 
3(3)  ,PROB<  15)  ,  XNC650),  YN<65u  ),  2N<650 ),  Th<650),  TXC650) 

4 , T Y ( 650 > , TZ  <  650 ) 

EQUIVALENCE  <KM.US<21 1 >), <KM,MK), <NP,US<201 )), <SVM,WS<210)) 
1,CV2,Y3),<n!X,U$<205)),<MIV,US(206)),<MlZ,WS<2C7)) 

INTEGER  PAGE 
INTEGER  SVH 

5  FORMAT (49H0XVZ  POTENTIAL  FLOU  PROGRAM  SECTION  5,  VERSION  4  ) 

WRITE  <6, 5) 

C  ft.  READ  INPUT 

READ (5, 25)  NOBP, I  EDITS, I  READ 
C  A.  READ  THE  OFF  BODY  POINTS 

NOB=NOSP 
DO  10  1=1, NOB 

READ <5, 26)  XP< I >,YP<I >,  ZPC I ) 

IF  (EOF <5 >  .EC!.  0.  )  GO  TO  10 
NQBF- 1 - 1 

WRITE <6, 9)  NOBP, NOB 

9  FORMAT ( 1 HO, I5,31H  OFF  BODY  POINTS  SPECIFIED  NOT  ,13) 

GO  TO  11 

10  CONTINUE 

11  CONTINUE 

25  FORMAT <314 ) 

26  FORMAT <3F  12. 5) 

P=1. 

READ  <03>  (PFlOBd  ),  1  =  1,  15) 

WRITE  (6,90)  <PROB<!  ), 1  =  1, 15) 

90  FORMAT < 1H0, 18A4) 

WR I TE<6, 9 1 )  NOBP, I  ED  ITS, I  READ 

91  FORMAT (SHONOBP  =',14  /8H  IEDIT5=,  14/SH  I  READ  =,14) 

WR I  TEC 6,92) 

92  FORMAT <  17H0  OFF  BODY  POINTS  /  4H  PT. ,  5 IX, 1HX,  12X,  1HV, 12X, 1H2) 

WR  I TEC6, 93 )  <I,XF'<I  ),VP<I  ),ZF'(I  ),  1  =  1, NOBP) 

93  FORMAT C IX .  1 3 , 2X, 3F 13 . 5  > 

C  B.  READ  THE  PARAMETERS,  T  ARRAY  AND  SOURCE  FROM  TAPE  31 

READ  <03 )  <US<I ), 1=1,220) 

READ  (03)  <TX< I >,TV(I ),T2<! ),XN(I ),VN<1  ),2N<I ),TACI ), l  =  1,NP) 

READ  <03 )  SKIP 
C 

C  FORMERLY:  WS<220)  ,E0.  2. 

C 

I F  < WS  < 220 )  . EQ .  5. )  READ <  03  >  SK I P 
READ  <03 )  (SKI  ),  l  =  1,NP) 

READ  <03 )  <S2< I  ), 1  =  1, NP ) 

READ  (03)  <S3<  I  ),  l  =  1,NF  ) 

C  C.  READ  THE  FIRST  BLOCK  OF  THE  B  ARRAY 

READ  (04)  <B< I  ), 1  =  1,241) 

K=  1 
J=1 

DO  100  1=1, NOBP 

C  D.  SET  THE  PARTIAL  VELOCITY  TO  THE  FREE  STREAM  VELOCITY 
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UXKi  >=-1.0 
WKI  >=0. 

UZKI  >=0. 

UX2< I  )=0 
UV2< I >=-10 
UX3CI >=0.0 
UV3< I  >»0 . 0 
UZ?< I >=-1.0 

100  UZ2( I )=0 . 

290  IFCBC 1 >-P >291,295, 291 

291  WRITE  (6,98)  B(1>,P 

98  FORMAT <28HG  POINTS  OUT  OF  ORDER  B<  1  >=, 1F4  0,4H  P=,1F4.0> 
STOP 

E  START  LOOP  OUEF;  THE  QUADS. 

295  J=2 

FI  PICK  UP  QUAD.  INFORMATION 

296  X1=B(J> 

V1=B(J+ )> 

X2=B(J+2> 

V2=B< J+3 > 

X3=B< J+4 > 

X4=B<J+5 > 

V4=B< J+6 > 

XC-TX(K > 

VC=TV(K> 

ZC=TZ(K> 

A  =TA< K  > 

XX=B< J+7 > 

YX=B< J+3 > 

ZX=E:(  J+9  > 

XV=B<J+10> 

VY=E:(J+ 1 1  > 

ZV-BC J+ 12 > 

F2  COMPUTE  LENGTH  OF  SIDES  OF  QUAD. 

D 1 2=SQ2F ( X 1 , X2 , V 1 , V2 , 0 . , 0 . > 

D23=SQ2F<X2,X3,V2,V3, .0, ,0> 

D34=S02F < X3 , X4 . V3 . Y4 , . 0 , . 0  > 

D4 1  =S02F  <  X4 . X 1 , V4 . Y 1 , .0. .0) 

F3  COMPUTE  SLOPE  OF  SIDES 
I F  X2-X3  >305 , 300 , 305 

300  Cl 23=1. 

GO  TO  310 

305  CM23=  C  V2-'v  3  >  /<  X2-X3  > 

Cl  23=0. 

310  I F1X3-X4 >3 15, 311,315 

311  C 1 34= 1 . 

GO  TO  320 

3 1 5  0134= < Y4-V3 > / < X4-X3  > 

C i 34=0 . 

320  IFCX4-X1  >325,321,325 

321  C 1 4 1=1. 

GO  TO  330 

325  CM4i= r  v j -V4 >  f / x i -X4  > 

C 1 4 1=0 

330  I F  O*  1— >335 , 33 1 , 335 

331  C 112=1. 

GO  TO  340 

335  CU12=<Y2-Y1 >/(X2-X1  > 

Cl  12=0 

F4  COMPUTE  OUAORAPOLE  MOMENTS 

340  CIXX=B(J+ 17> 

CIXV=B(U+ 18 > 


C!VY=B(J+1y) 

F5  COMPUTE  SIN  AND  COS  OF  SLOPE  ANGLE  FOR  EACH  SIDE 
CY12=<Y2-Y!)/D12 
CV23=  (  V3-Y2  )  /D23 
CY34=(Y4-Y3VD34 
CY4 1=(Y 1-V4  )/D4 1 
CX 12=(X 1-X2  >/D 12 
CX23=  ( X2-X3 ) /D23 
CX34=CX3-X4 )/D34 
CX4 WX4-X 1 >/D4 1 

F6  COMPUTE  MAX  D ! AG ! UAL 
ST =SQ2F < X 1 , X3, V 1 , V3, 0  . , 0 .  ) 

ST2-S02F (X2 . X4 , Y2 . Y4 . . 0 . .  CD 
I F ( ST-ST2  >34 1 . 342 . 342 

341  ST=ST2 

G  START  LOOP  OUER  THE  OFF  BODY  POINTS 

342  DO  530  JQ=1.N0BP 
I  S=  1 

XCQ=XP< JO > 

YCO=VFCJQ) 

ZCO=ZP(JO> 

J1=  1 

345  RPQ=SQ2F  <  XC ,  XCQ ,  VC ,  VCQ ,  2C  SCO  > 

H.  DETERMIN  METHOD 

I F  < RPQ-ST*4  >350 ,350,460 
350  X=< XCQ-XC  >*XX+  CYCC'-YO )*YX+  ( 2C0-2C )+2X 
V= (XCQ-XC >*XY+( YCQ-VC >*VY+(ZCQ-ZC >+ZV 
2= < XCQ-XC >*XN(K )+ (YCQ-VC >*VN(K )+<2CQ-2C  >*Z?KK ) 
IF(RPQ-ST*2 . 5 )355 . 355 , 400 

I,  COMPUTE  I  HOUSED  UELOCITY  BY  EXACT  METHOD 
355  R1»SQ2F(X,X1,V,V1,Z.O. ) 

R2=S02F(X,X2JY,V2JZ,0. > 

R3=SQ2F(x!x3.y'v3.z'o. > 

R4=S02F  < X ,  X4 ,  Y,  Y4 ,  S , 0 . > 

IF((R1+R2).LE.D12>  GO  TO  1000 
I F< (R3+R2 ) . LE . D23 )  GO  TO  1000 
IF((R3+R4>.LE.D34>  GO  TO  1000 
IFC <P 1+R4 >.LE. 04 1 >  GO  TO  1000 
CLR 1 =RLOG( <R 1+R2-D 12  > / < R 1 +R2+D 12  » 

CLR2=RL0G( (R2+R3-D23 >/(R2+R3+D23  )  > 

CLR3sRwOGc (R3+R4-D34 )/(R3+R4+D34 ) ) 

ClA4=PL0G< (R4+R 1-D4 1 >/<R4+R 1+D4 1 ) > 

TUX=CV 12*CLA 1+CY2?*CLfl2+CY34*CLfi3+CY4 1*CLR4 
TUV=CX 12+CLfi 1+CX23*CLR2+CX34*CLR3+CX4 1*CLR4 
TUZ=0 , 

i F ( ABS <Z>-. 00 1*ST >375 ,361. 35 1 
361  ZSQ=2**2 

E1=ZS0+(X-X1 >++2 
E2=2SQ+(X-X2  y**2 
E3=2SQ+ ( X-X3  y**2 
E4=2SQ+  < X-X4 )*+2 
H 1=< V-V 1 )*<X-X  1 ) 

H2=(Y-Y2>*(X-X2) 

H3=  <  Y- V3 )* ( X-X3  > 

H4=(V-Y4 >*(X-X4> 

IF (Cl  12)363,353.354 
363  US  1  =  (Cr1 12+E  1-H  t  >/(Z+R  1  > 

MS  2=  <  CM  12+ E2-H2  >  /  ( Z+P.2  > 

AT 1=ATAM(US  1 ) 

AT2=ATAN(US2 ) 


TUZ=AT  1-AT2 

364  I F(C 1 23 >366, 366, 367 


366  RT3=ATRN':  CCM23+E2-H2  >/<Z*R2  )  > 

RT4=ATRH< <CM23*E3-H3 >/<Z*R3 >  > 

T0Z=TUZ+RT3-flT4 

367  I FCCI 34 >368, 368, 369 

368  RT6=RTFlN  ( <  CM34*E3-H3  >  /  < 2*R3  >  > 

AT6=RTRfi<  (CM34*E4-H4  >/<Z*R4  > ) 

TUZ=TUZ+AT5-RT  6 

369  I FCCI41 >370,370,375 

370  RT7=fiTRN ( ( CM4 1 +E4-H4  > } i 2*F:4 ) ) 

RT8=RTHN( <CM4 1*E 1-H 1 >/<Z+R 1 > > 

T0Z=T0Z+AT?-RT8 

375  GO  TO  450 

C  J.  COMPUTE  I  MOUSED  VELOCITY  BY  OURDRRPOLE  METHOD 

400  RF’Q3=F;PQ*+3 

RPQ7=  ( RPQ3**2  >*RPO 
US1=  X/RFQ3 
XS0=X**2 
V$0=V**2 
ZS0=Z**2 

PS=V£0+ZSu-4 . +XSQ 
0$=XSQ+2S0-4  ♦YSO 
US2=X*  ^9 .  *PS+3Q .  *XSQ  >/F;FQ7 
WS3=3 .  *V+Pc>/FiF'Q? 

US4=3 *X*QS/RPQ7 

TUX=A*US 1-C I XY+0S3-C I XX*US2-C I YY*U£4 
WS1-V/RPQ3 

W$2=V* < 9 . +QS+30  *YSQ > /RP07 
TUV=h  •  „S  *. -3 ;  XX^SS-C ;  XV*yS4-C  i  VV*»S2 
TUZ=Z  *  <  ft  /'FiF'03-3 .  *<C  I  XX*P$-5 .  *C !  XV+X+V+C  I  VV*QS  >/RPQ7 
450  UX< ! S  >=T0X*XX+TUV*XY+TU2*XN ( K > 

W<  I S  >=TOX* VX+TUY*YV+7UZ*YN<K  > 

uzc :  s  >=tvx+zx+toy*zy+toz*znck  > 

GO  TO  470 

C  K.  COMPUTE  I  HOUSED  UELOCITV  BY  MONOPOLE  METHOD 

46 C  RF!PQ3=fl/CRP0**3> 

t.JX  (I  S  >=  <  XCO-XC  >*RRF’Q3 

UY  ( I S  )=  (VCO-VC  >+RRF'03 
UZ( !  S  >=  (ZCQ-ZC  >*RRPQ3 

C  L  REFLECT  OFF  BODY  POINT  IN  PLANE  OF  SVMETRV 

470  00  TOf 4S0, 485, 490, 495, 500,505, 5 10 . 5 15 >, IS 

480  UKJ1  )=0X(  1  > 

U2(J1)=UXa> 

03  <  J 1  )=U\Y  1  > 

0  KJ1+  1  >=UV(  1  > 

02(01+  1  >=UY( 1 > 

US (J1+ 1  )=UV< 1 > 

OKU  1+2 >*0Z<  t  > 

U2(J1+2)=0Z( 1 > 

U3‘.  J 1+2  >=IJZ(  1  > 

I F < SVM >  525,525,481 

481  1 8=2 

C  XZ  SVMETRV 

VCQ=-7C0 
GO  TO  345 

485  I F ( SVM- 1 >5 1 7 , 5 1 7 , 486 
0  XV  SVMETRV 

4S6  IS-3 

zco—zco 

GO  TO  345 
490  1 8=4 

VCQ=-YCO 
GO  TO  345 
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495  IFcSVM-2  )5 16, 5 16,436 

YZ  SYMETRY 

496  !S=5 
XCQ=-XCQ 
GO  TO  345 

500  IS=6 
YC0=-YC0 
GO  TO  345 
505  !£=? 

2C0=-ZCQ 
GO  TO  345 
510  I S=8 

VC0=-VC0 
GO  TO  345 

M.  ADD  COHTP 1  BUT  I ONS  OF  ALL  REFLECT! OHS 
515  0  K  0 1  >=0  K  J  1 >+0X < 8  >+0X ( 7  )+UX ( 6  >+UX <  5  > 

U2< J 1 )=U2< J 1 >-0X(8>+0X<? >+UX(6  >-0X<5  > 

03  <:  J 1  >=U3  <  J 1  >-0X  ( 8  >-UX  <  7  >+UX  ( 6  >+0X  <  5  > 

0 1  ( J 1+  1  )=U  K  J 1+ 1  )-UV  <  8  >+UV  ( 7  >+UV  ( 6  )-UV(5  > 
02 J 1+ 1  >=02(0 1+  1  >+0\XS  >+7V(7  >+NY<6  >+UV<5 ) 
U3 < J 1  + 1  )=03  ( J1+  1  >+ny  <  0  -ny  r  7  >+07  < 6  )-UY < 5  ;> 
UKJ1+2  >=0  KJ1+2  >-UZ  <  8  >-UZ  ( 7  >+02  <  6  )+UZ  ( 5  ) 
02  < J 1+2  >=02< J 1+2  >+0Z< 8  >-02  < 7  >+UZ<6 >-0Z<5 > 
03(01+2 >=03 ( J 1 +2 >+UZ<8 >+0Z<?  >+0Z(6 >+0Z(5 > 
515  01 <  J 1  >=U  H  J 1  >+UX  <  4  >+UX  <  3  > 

02  <  J 1  >=02  <  .J  1  >+0X  (  4  >-0X  <  3  > 

03 < J 1 >=03 < J 1 >-UX < 4  >-0X < 3  > 

0 1 C J 1+ 1  )=U  1  < J1+ 1  )+UV(4 >-UV<3 > 

02 <01+1 >=02 <01+1 >+0V <  4  >+0V <  3  > 

03 ( 01+1 >=03 <01+1 >-0V <  4  >+0V < 3  > 

0 1 <01+2  >«U 1 <0 1+2  >-0Z<4 >-0Z<3> 

<'2  <01+2  >=02  <  0 1  +2  >-0Z  <  4  >+02  <  3  > 

03 <0 1+2  >=03 < J 1+2  >+0Z<4  >+0Z<3 > 

517  01  <01  >=01<01 >+0X<2 > 

02  < 0 1  >=02 < 0 1 >-0X  <  2  > 

03(01 >=03(01 >+0X< 2 > 

01 <01+1  >=01 (01+1  >-07(2) 

02 <  01+1 >=02 <  01+1 >+0V ( 2  > 

03  < 01+1 >=03 < 01+1 >-07  <  2  > 

0 1(0 1+2  >=0 1 (0 1+2  >+0Z<2  > 

02(0 1+2  >=U2(0 1+2  >-02 <2 > 

03 < 0 1 +2  >=03 < 01+2  >+0Z  < 2  > 

c-“-c  > 

OX  1 <00  >=0X 1 <0Q >+0 1 < 1  )*S 1<L  > 

07 1 <  00  >=0V 1 <  00 )+0 1 <  2  >+S 1 <  L  > 

OZ 1 ( 00 >=02 1 <  00  >+U 1 < 3  >*S 1 < L  > 

0X2(00  >=0X2(00  >+02 < 1 >+S2(L  > 

072  ( JQ  >=i  .'72 <  00  >+02  ( 2  >+S2  <  L  > 

022 < JQ  >=0Z2 < JO  >+02 ( 3  >*S2 < L  > 

0X3(00 >=0X3 <00 >+03 ( 1 >*$3<L > 

073(00 >=073(00 >+03<2  >*$3<L > 

023(00 >=023(00 >+03(3  >+53fL  > 

530  C0NT I NIJE 

N.  END  OF  LOOP  OOER  OFF  BODY  POINTS 
585  P=P+ 1 

K=K+ 1 
J=J+20 

) F  < K-NF  >586 , 586 , 599 
5c:6  i  F  <  0—24 1  >296 , 590 , 5 : 

0  READ  NEXT  BLOCK.  OF  B  ARRAY  IF  NEEDED 
590  READ  <04  >  <B< ' >, 1  =  1 . 24 1 > 


IF(B(1)-P  >291.296,231 

P.  END  OF  LOOP  OUER  QUADS 
599  CONTINUE 

PAGE  =  1 

IF  < IEDIT5  .EQ.  1>  GO  TO  825 

601  FORMAT <4H  PT. ,  1 1X, 1HX, 12X, 1HV . 12X. 1HZ,  14X,2HUX,  1 1X,2HUV, 1 1X,2HUZ 
1,  14X,2HCP> 

602  FORMAT (?H  X  FLOW) 

603  FORMAT C7H  V  FLOU) 

604  FORMAT <7H  Z  FLOW) 

605  FORMAT (1H1, 15A4, 10X, 15H0FF  BODV  POINTS  ,  10X.5HPAGE  , IS) 

IF  (MIX.EQ.O)  GO  TO  700 

WR I  TEC 6. 605)  PROB,PAGE 
WRITE  <6,602) 

WRITE  <6,601) 

L I NE= 1 
LAST=53 

606  I F (NQBP-LAST  )607 ,610,610 

607  LAST=NOBP 

610  DO  615  l=LINE,LAST 

6 1 1  FORMAT < 1X.1I3.2X. 3F 13 . 5 , 2X . 3F 13 .5.3X.F 13 . 5 ) 

Q.  COMPUTE  PRESSURE  AND  EDIT  3  BASIC  FLOWS 
CP  1= 1 . -<UX 1 < I >**2+UV 1 < I >**2+UZ 1 < I >**2 > 

615  WRITE  (6,611)  I  ,XP<I  ),VP<I  ),ZP<I  ),UX1<!  ) ,  UV 1  <  I  ),UZH!  ).CP1 
L INE=LAST+ 1 
LRST=L INE+54 
PAGE=FAGE+ 1 

I F  <  L I NE-NOBP >620 , 620 , 700 
520  WR ! TE  < 6 , 505 >  PROS , PAGE 
WRITE  <6,601) 

GO  TO  606 

7  op  |p  <M!V  EM  n>  RM  jn  AMR 
WRITE <6 ,605 )  PROS , PAGE 
WRITE  (6,603) 

WRITE  (6,601) 

L I NE=  1 
LAST =55 

706  I F  <  NOBP-LAST >707 ,710,7 10 

707  LAST=NOBR 

710  DO  715  I =L I NE , LAST 

CF2=1 . -<UX2< I  )**2+UV2< I >**2+022 < 1 >**2> 

7 15  WR i TE  (6,611)  I , XP< I > , VP< I ).2P(I )  UX2( I  ) , UV2< I > , U22< I >, CF2 
L INE=LASf+ ! 

LAST =L I NE+54 
PAOE=pmOE+ 1 

I F  <  L I NE-NOBP >720 . 720 , 800 
720  WRITE <6, 605)  PROS, PAGE 
WRITE  (6,601) 

GO  TO  706 

800  IF  (MIZ.EO  0)  00  TO  825 
WRITE (6, 605)  PROS, PAGE 
WRITE  (6,604) 

WRITE  (6,601) 

L INE=  1 
LAST =55 

805  I F (NOBP-LA: T >807, 810, 8 10 

S07  LAST=NOBP 

810  DO  815  I =L I  HE, LAST 

CP3=  1  -<UX3>:  I  >**2+UV3<  I  >**2+UZ3(  I  >**2  ) 

S 15  WR I TE  (6,611)  I , XP < I  ) , VP < I > . ZP < I  ) , UX3 1  I  ) , UV3( I  ) , UZ3< I ) , CPS 
L I NE=LRST+ 1 
LAST =L I NE+54 


PAGE=PAGE+ 1 

I F ( L I NE-NOBP >820 , 820 , 825 
i  UR I  TEC 6, 605)  PROS, PAGE 
URITE  (6,601) 

GO  TO  806 
i  J  =  1 

.  IF  (IRERD.EQ  0)  GO  TO  827 
READ (5, 26)  UX4,UV4,UZ4 
IF  < EOF (5 ) . NE  0.)  GO  TO  QOO 
GO  TO  828 
'  UX4=US(J) 

UV4  =  USCJ+20) 

UZ4  =  US ( J+4G ) 

:  CP=UX4*+2+UY4**2+UZ4**2 
I F< CP >900, 900 ,830 

R.  COMPUTE  FOURTH  FLOU  AMD  EDIT  IT 
i  L I NE= 1 
LAST=5 1 

URITE (6, 605)  PFlOB , PAGE 

FORMAT ( 10HGCNSET  FLOU  UX  =  ,F7.3/15X,4HW  =, 
1F7.3/15X.4HU2  =,F7.3> 

URITE  (6,831 )  UX4 , UV4 , UZ4 
URITE  (6,601 ) 

5  I F (NOBP-LRST )837, 840, 34G 
’  LAST=NO£P 
)  DO  845  l=LINE,LRST 
U.XF--U.X4+UX  1  ( I  )-UV4*UX2< !  >-UZ4*UX3<  I  > 

UVF--UX4+UV 1  ( I  )-UV4*UV2( I >-UZ4*UV3< I ) 
UZF--UX4+UZK  I  >-UV4*UZ2C  I  >-UZ4*UZ3<  I  ) 

CP4=  1 .  - <  UXP**2H'YP**2+UZP**2 ) /CP 
I  URITE  (6,611)  I  , XP ( I  ) , VP ( I  ) , ZP <  I  > , UXP , UVP , UZP , CF'4 
L !NE=LRST+ 1 
LftS7=L INE+54 
PR0E=PAGE+ 1 

I F ( LI NE-NOBP )850 , 850 , 860 
i  URITE (6, 605)  F'ROE , F'RGE 
WRITE  (6,601) 

GO  TO  835 

i  J  =  J+1 


i  URITE(6, 1001)  JQ,L,XP(JQ),VP(JQ),ZP(JQ) 
i  FORMAT ( 16H00FF  BODV  POINT  , I3,23H  ON  BOUNDARV  OF  OUAD 
1  3H  X=,F 12. 5,5X, 2HV=,F  12 , 5,5X,2HZ=,F 12.5) 

GO  TO  530 
i  CONTINUE 

S  FlEUIND  TAPES  AND  STOP 
REUIND  03 
REWIND  04 
STOP  5 
END 

FUNCT I  ON  SQ2F ( X 1 , X2 , V 1 , V2 , Z 1 , Z2 ) 

X=X  1-X2 
V=V  1-V2 
Z=Z  1-Z2 

RS=Z**  2+ Y**  2+  X*+2 

Fi=RES  ( X  )+RBS  ( V )+ ABS  ( Z )+  1  0E-20 

R=Fl+RS/R 

R=  25+fi+RS/R 

R=  R+RS/R 

SQ2F=  .  25*R+P.S/R 

RETURN 

END 


o  o  o  o  o 


APPENDIX  U1  -  XVZPF  SECTION  PF6 


PROGRAM  PFP6C I NPUT= 128, TAPE  16, GUTPUT= 128, TAPE03, TRPE04, 

1  TflP£5= I NPUT, TAPE6=0UTPUT, TAPE3=TAPEG3, TAPE4=TAPEG4 ) 


XVZ  POTENTIAL  FLOW  PROGRAM  VERSION  4  SECTION  6 
COMPUTES  VELOCITIES  AND  PRESSURE  COEFFICIENTS  FOR 
OFF  BODV  STREAMLINES 


COMMON  XP<  100),YP(  100>,ZP<  100>,VX1(10G),US<220> 

1 ,  UV 1  ( 1 00  > ,  VZ  K  1 00  > ,  UX2  ( 1 00  > ,  UV2  <  1 00  ) ,  UZ2  (  100),  UX3  (  1 00  ) ,  UV3  (100) 

2, VZ3<  100),VX(8),  VY(8),  VZ(8),S1(650),S2(650),S3(650),  V1(3),  V2(3),V3 
3(3 )  , PROB( 15 )  , XN(650 ) , S'N<650 > , ZNC65G ) . TX(65C ), TV(650 ) . TZ(650 ) 

1,B(  13000),  XT  ( 100>,YT(  100),ZT(  10b>,AP(5).GM(4),SKV(  100), 

1  SKZ ( 100),SKX( 100), DX( 100), DV( 100), DZ( 100), CP( 100)  ,TA(650) 

EQU I  VALENCE  (KM , WS(2 11 ) ), (KM , UK ), (NP, US (20 1 ) ) , CSV”’  WS(2 10 ) ) 

1 , (Y2, Y3 ) 

INTEGER  SYM  ,P 

A.  READ  INPUT 
WRITE  (6,5) 

6  FORMAT (3 1 4, 4F 12. 5) 

8  FORMAT (3F 12.5) 

5  FORMAT (4QH0XYZ  POTENTIAL  FLOW  PROGRAM  SECTION  6, 

7  FORMAT  (1X.9F12.5) 

20  FORMAT C  1H1, 14.31H  STREAMLINES  TO  BE  COMPUTED  AT 

1  ,F8.4,28H  T  FOR  AN  ONSET  VELOCITY  OF  ,3F8.4) 

2 1  FORMAT (IX, 1 5HST ART  I NG  Pu I  MTS / , 3X , 2HPT ,5X, 1HX, 1 1X, 1HV, 1 1X, 1HZ) 

22  FORMAT ( IX, 14, 3F 12 . 5 ) 

READ  (03)  (PROB(! >, 1=1, 15) 

WRITE  (6,90)  (PROS (I ), 1=1, 15) 

90  FORMAT ( 1H0, 15A4) 

B.  READ  THE  PARAMETERS,  T  ARRAY  PND  SOURCE  FROM  TAPE  31 
READ  (03)  (US( I ), 1=1,220) 

READ (03 >  (TX< !  >,  TV C I ),  T2<  /  >, XN<  I  >,  YH< !  >, ZN(  l  >,  Tft(  I  >,  1  =  1, NP) 

READ  (03)  SKIP 

READ (03)  SKIP 


VERS  I  ON  4  ) 

,14,  1C.H  STEPS  OF 


IF  (  US (220)  ,E0.  2.  ) 

READ  (03)  (SKI  >,  1  =  1, NP) 
READ  (03)  (S2( I  ), 1  =  1, HP) 
READ  (03)  (S3(l ), 1=1. NP) 

C.  READ  THE  E  ARRAY 
UZ  =  NP 

WZ=(WZ+ 11.0)/ 12.0 
NB  =  UZ 


IF=24 1 

DO  12  IP  =  f,NB 

READ  (04)  P,  (B(  I ),  1  =  1  S.,  IF) 

IS= IS+240 
12  I F= IF +240 
AP< 1 >=  .5 
AP(2)  =  .5 
AP ( 3 )  =  1 . 

AP  ( 4  )=0 . 

APC5)=0 
GM( 1 )  *  1./6. 

GM(2)  =  1./3. 

GM(3)  =  1./3. 

83  READ  (5,6)  NOBP , NST , ! END , DT , VX I , VY I , VZ 1 
VSQ=VX I ++2+VV I ++2+VZ I **2 
HGB=NGBP 
DO  10  1=1, NOB 

READ(5,  8)  XP(I ),VP(I ),ZP(I ) 


i 

■ 


>■*  V^^'vTV^VW  'J*  '.’g’.^'TVVVl’gVWI  l'f  m  i.-M  yr  f  IT  l-ty*  r»  t 


'.’’V'VTl''  rr*  L’’  rnmir*  v-»  r-> 


IF  (E0FC5)  .EQ.  0. )  GO  TO  10 

N0BP= I - 1 

WR I TE ( 6 , 9 )  NOBP, NOB 

9  FORMAT < 1H0, 15,  28H  STREAMLINES  SPECIFIED  NOT  , 13) 

GO  TO  11 

10  CONTINUE 

11  CONTINUE 

WRITE( 16)  N0BP,NST, IEND,UXI , Wl ,UZI 
UR  I TE(6,20 )  NOBP, NST, DT, UX I , UV I , U2 1 
UP.  t  TE<6, 2 1  > 

I4RITE<6,22)  <  I ,  XP<  I  ),  VP<  I  >,2P<  I  ),  1  =  1, NOBP) 

C  NGBP  -  NUMBER  OF  STREAMLINES  TO  BE  TRACED. 

C  NST  -  NUMBER  OF  STATIONS  AT  UHICH  STREAMLINES  SHOULD  BE  COMPUTED. 
DO  15  1=1, NOBP 
XT < I  )  =  XP(I ) 

VT< I >  =  VP< I  ) 

ZT< I )  =  ZPC I  ) 

SKXCI )=0 . 

SKV< I )  =  0. 

15  SKZCI  )  =  0. 

I  TC=0 
I  RK=5 
98  K=1 
P  =  1 
J=1 

DO  100  1=1, NOBP 

C  D.  SET  THE  PARTIAL  UELOCITV  TO  THE  FREE  STREAM  UELOCITV 

UX  1  <  I  )=- 1 . 0 
UV  1  <  I  )=0 . 

UZKI  )=0 . 

UX2< I  )=0 . 

UV2< I  )=- 1 . 0 
UX3CI  )=0.0 
UV3( I )=0. 0 
UZ3< I  )=- 1 . 0 
100  UZ2< I )=0. 

C  E.  START  LOOP  OUEFl  THE  QUADS. 

295  J=2 

C  FI  PICK  UP  QUAD.  INFORMATION 

296  X  1=B(J) 

V  1=B( J+ 1  ) 

X2=B(U+2> 

V2=B ( J+3 > 

X3=B< J+4 ) 

X4=E'(  J+5  > 

V4=B(U+6 ) 

XC=TX(K ) 

VC=TV<K > 

ZC=T2(K ) 

A  =TA(K) 

XX=B( J+7 ) 

VX=B< J+8 > 

ZX=B( J+9  > 

XV=B  < J+ 1 0 > 

VV=B< J+  1 1 ) 

ZV=B< J+ 12) 

C  F2  COMPUTE  LENGTH  OF  SIDES  OF  QUAD. 

D 1 2=SQ2F  O'  1 , X2 , Y 1 , V2 , 0 . , 0  > 

D23=S02F < X2 , X3 , V2 , V3 , . 0 , . 0 ) 

034=SQ2F  CX3 , X4 , V3 , V4 , . 0 , . 0  > 

04  1  =S02F ( X.4 , X 1 , V4 , V 1 ,  .0.  .0) 

C  F3  COMPUTE  SLOPE  OF  SIDES 


1  17 


■  ^  -  a  - T  i 
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I F  F  X2-X3  >305 , 300 , 305 

300  C 1 23= 1 
GO  TO  310 

305  01123=  <  V2-V3 )  /  <  X2-X3 ) 

Cl  23=0. 

310  IF(X3-X4 >315,31 1,315 

311  Cl  34-1. 

GO  TO  320 

315  CM34=<Y4-Y3  >/FX4-X3 ) 

Cl  34=0 

320  I FFX4-X  1)325, 32 1,325 

321  C 1 4 1=  1 . 

GO  TO  330 

325  CM4 1=<V 1— V4 >/<X 1-X4  > 

Cl  4 1=0 

330  I F ( X 1 -X2  >335 ,331, 335 

331  Cl  12=1. 

GO  TO  340 

335  CM 12=< Y2-Y 1 >/(X2-X 1 > 

Cl  12=0. 

F4  COMPUTE  QURORRPOLE  MOMENTS 

340  C I XX=B<J+ 17  > 

C I XV=B< J+ 18) 

C I YV=B ( J+ 1 S  > 

F5  COMPUTE  SIN  FIND  COS  OF  SLOPE  fiNGLE  FOR  EACH  SIDE 
CV  12=<V2-V  1  >/D  12 
CV23=  F V3-V2  >/'D23 
C V34=  <  Y4-V3  > /034 
CV4 1  =  F  V 1  -  V4  >  /D4  1 
CX 12= (X 1-X2  >/D  12 
0X23=  < X2-X3  > /D23 
CX34= (X3-X4  >/034 
CX4 1=FX4-X 1 >/D4 1 
F6  COMPUTE  MAX  0  IPG  INAL 
ST =S02F<X  1 , X3 , V 1 , V3 , 0 . ,0. > 

ST2=SQ2F ( X2 , X4 , V2.V4,  0,0) 

I FFST-ST2 >34 1 , 342  342 

341  ST=ST2 

0.  STAFtT  LOOP  OvEFi  THE  OFF  BODY  POINTS 

342  DO  530  JQ=1,NQBP 
IS=  1 

XCO=XP< JQ > 

VCQ=V=‘,JO> 

ZUU=Zr  FJu ) 

J1=1 

345  RPQ=SQ2F <XC . XCQ . VC . VCQ, ZC, ZCQ ) 

H.  DETERM  IN  METHOD 
I F<RPQ-ST*4 >350, 350 . 460 

350  X=FXCCi-XC )*XX+< VCQ-YC  >*YX+<ZCQ-ZC >*2X 
V=FXCQ-XC  >*XY+F VCO-VC >*VY+FZCQ-ZC )*ZV 
Z=<XCO-XC  >XN(K  >+CVCG-VC  X+VNFK  >+<ZCQ-ZC  >*ZN<K> 

I F<RP0-ST*2 . 5  >355,355,400 

I  COMPUTE  I MOUSED  VELOCITY  BY  EXACT  METHOD 

355  F; !  =SQ2-  •.  X ,  X 1 ,  V ,  V 1 ,  Z ,  0  > 

R2=302F  <  X , X2 , V , Y2 , Z , 0 .  > 

R3=SC2F ( X , X3 , V , V3, Z , 0 .  > 

R4=S0'2F CX,  X4  ,  Y  ,  Y4  ,2,0.  ) 

CLP  1= ALOOF (R1+R2-D 12  >/<R 1+R2+D 12  >> 

CLR2=AL00F FR2+R3-D23 )/<R2+R3+D23 > > 

CL A3= ALOOF  FfiU+R4~D34  )/FRo'+R4+D34  >  > 

CLR4= ALOOF  FR4+R1-D4 1  >/FP,4+R  1+04 1  > ) 

TVX=CV  12:fCLA  »+CV23*CLft2+CY34*CLA3+CV4 1*CLP4 


TUV=CX  1 2*CLfl 1 +CX23*CLR2+CX34*CLR3+CX4 t*CLR4 
TU2«0 . 

1 F (RBSCZ )- . 00 1*ST >375, 36 1 , 36 1 
230=2**2 

E  1=ZSQ+<X-X 1  >**2 
E2=ZSG+<X-X2>**2 
E3-2SQ-KX-X3 >**2 
E4=ZSQ+(X-X4 >**2 

Hi=<v-vi>*a-xi> 

H2=(V-V2 >*<X-X2> 

H3=(V-V3 >*<X-X3> 

H4= ( V-V4 )*  (  X-X4 ) 

IF(C1 12 >363, 363 ,364 
US1=(CM12*E1-H1 >/(Z*R1 > 

WS2=(CM 12*E2-H2 >/<2*R2 > 

RT1=R7R?KWS1> 

RT2=flTRN<US2 > 

TUZ=RT 1-RT2 

!F (Cl  23 >366, 366, 367 

RT3=RTRN< <CM23*E2-H2  >/<Z*R2  >  > 

RT4=RTfiM<  <CM23*E3-H3  >/<Z*F;3  >  > 

TUZ=TUZ+fiT:?-ftT4 
IF (Cl  34 >363, 368, 369 
RT5=flTRN(  (CM34*E3-H3  >/<Z*R3  >  > 

R75=RTRN(CCM34*E4-H4  >/(2*R4  >  > 

TU2=TU2+RT5-RT6 
l  C(C! 41 >370, 370, 375 
RT7=RTRH(  (014 1*E4-H4  )/<Z*R4  >  > 
ftT8*RTftm  <Cr*4  >/<2*R  1  >  > 

T\IZ=TUZ+RT7-RT8 
GO  TO  450 

J  COMPUTE  !  MOUSED  UELOCITV  BV  OURDRRPOLE  METHOD 
RPQ3=R?Q**3 
PPQ?-(P,PQ3**2  >*R PQ 
US  1 =  X/RPQ3 
XSG=X**2 
S'S0=V**2 
ZSQ=Z**2 

F3=VbO+ZSQ-4 .  *XS0 
CS=XSCi+ZSQ-4 .  *VSC 
US2=X*C5 .  *F'S+30 .  *XSQ  )/RPO? 
ljF3*3  *m*ps  /pPQ7 
US4=3  *X*(.ib  /RPQ7 

TUX=fi*U3 1-C I  X7*US3-C  I  XX*US2-C  I  VS'*U$4 
US  1=V/RP03 

US2=V*C9 . *GS+3G . *VSQ  >/RPQ7 

TW=R*US  1-C I XX+US3-C !  XY*US4-C I VV+US2 

TU2=Z* ( fl /RPQ3-3 . * < C I XX*PS-5 . *C I XV*X*Y+C I W+QS >/RPG7 > 

UXC I S  >=TUX*XX+TUY*XY+TUZ*XN<:K  > 

UV<  i  S  > = T  U  X  *  V  X  ♦  T  U  V  *  V  V  ♦  T  U  Z  *  V  N  ( K  > 

UZ<  IS  )=TUX*ZX+TW*ZV+TUZ*ZrKK> 

GO  TO  470 

K.  COMPUTE  I  MOUSED  UELOCITV  BV  MOMC'POLE  METHOD 
RFlF'Q  :;sfi  /  t  RPu**y  > 

UX  ( I S  >=  <  XC'O-XC  >*RRP03 

l  IV ( ! s ' -  ( VCQ-VC >*RRP03 

UZ •.!£>=  (ZCQ-ZC  >*RRPQ2 

L  REFLECT  OFF  600 V  POIMT  IM  PLRflt  OF  SVHETRV 
GO  TO  < 480 , 485 ,450,495, 500 , 505 , 5 1 0 , 5 1 5  > ,  I S 

U1CJ1  >=UX( 1 j 
U2CJ1>=UX< 1 > 

U3f  J 1  >=UX',  1  > 


01(01+1  >=0V<  t ) 

02CJ1+1  >=W<  1  > 

V3CJ1+1  >=W<  1 ) 

U 1  (0 1+2  >=UZ(  1 ) 

U2C0t+2>=UZC1> 

U3(J1+2  >=0Z(  1 ) 

IF<SVM>  525,525,481 

4S1  IS=2 

X2  SVMETRV 

VCQ— VCQ 
GO  TO  345 

485  I FCSVM-1 >517,517,486 

XV  SVMETRV 

486  IS=3 
2CQ=-ZCQ 
GO  TO  345 

400  IS=4 
VCO=-VCQ 
GO  TO  345 

405  I F (SVM-2  >516,5 16,496 

VZ  SVMETRV 

4Go  I S=5 
XCC!=-XCQ 
GO  TO  345 

500  I S=5 
VCQ=-VCO 
GO  TO  345 

505  I S=7 

ZC0=-2C0 
GO  TO  345 

510  1 3=8 

VCQ=-VCQ 
GO  TO  345 

M.  ROD  COMTR I  BUT  I  OHS  OF  RLL  REFLECTION 

515  U  K  J 1  >=U K  J 1  >+UX ( 8 >+UX <  7  >+UX ( 6 >+0X ( 5  > 

U2<  J 1 >=02< J 1 >-0X<8  >+0X<7  >+UX<6  >-UX<5  > 

03 < j i >=03 C J 1 >-0X ( 8  >-UX  < 7  >+0X  < 6  >+0X < 5  > 

0 1(01+1 >=0 1(01+1 )-0V(8 >+OVC?  >+UV<5  >-UVC 
02  < 01+1  >=U2 <  J 1  +  1 >+0V  < 8  >+0V  < 7  >+0V  <  6  >+0V ( 

03  <  0 1 + 1 >=03  <  0 1 + 1 >+0V ( 8  >-0V ( 7  >+UV ( 5  >-0V  ( 

0 1  (0  1+2  >=0 1 (0  1+2  >-02 <8 )-0Z(7  >+0Z<6  >+0Z< 

02  ( J 1+2  >=02  ( J 1  +2  >+0Z  (8  )— UZC7  )+0Z(  6  >—02  ( 

03  C  J 1  +2  >=03  ( J 1  +2  >+UZ < 8  >+UZ ( 7  >+UZ ( 6  >+UZ ( 

516  0  K J 1 >=U U  J 1 >+0X C 4 )+0X  (  3 ) 

02  ( J 1  >=02  C J 1  >+0X  ( 4  >-0X (  3  > 

03 ( J 1 >=03  ( J 1 >-UX ( 4  >-UX ( 3 > 

0  K  J 1+ 1  >=0  KJ 1+ 1  >+0V<4  >-0V<3  > 

02(01+1  )=02  ( 01+1  >+0V  (.  4  >+0V  (  3  > 

03(01+1 >=03(0 1+ 1 >-UV(4  >+0V<3 > 

0 1  C0 1+2  >=0 1 (J 1+2  >-0Z(4 >-0ZC3 > 

02 C0 1+2 >=02 C0 1+2 >-UZ(4 >+0Z(3 > 

03C  J 1+2  >=03 (.01+2  >+0Z ( 4  >+UZ C 3  > 

517  01(01  >=01(01 >+0X< 2 > 

02 ( 0 1  >=02  < 0 1 >-UX ( 2  > 

03 ( 0 1  >=03 < 0 1 >+0X ( 2  > 

0 1  ( 01+1  >=0 '(01+1 >-07 c 2  > 

02(0 1+1  >=02(J 1+1>+0V(2) 

03CJ1+1 >=03  C 01+1 >-0V C  2  > 

01(0 1+2 >=01(0 1+2 >+0ZC2> 

02C0 1+2 >=02(0 1+2 >-UZ(2> 

03  C  0 1 +2  >=03C  J 1 +2  >+0Z ( 2  > 

525  OX  1  ■  00  >=0X  1  COCl  >+U  1  ( 1  >*S  1  (P  > 


in  in  in  in  in  in 


UY 1  (JQ  )=UY  K  JQ  )+U  K2  )*S  t  <P ) 

UZ  H  JQ  >=UZ KJQ  HU  K3  >*S  UP  > 

UX2CJQ )=UX2( JQ HU2( 1 >*S2(P ) 

W2  ( JQ  >=UV2  ( JQ  )+U2  <  2  >*S2  <  P ) 

U22< JQ >=UZ2< JQ HU2<3 >*S2(P> 

UX3  ( JQ  )=UX3  <JQHU3<  1  >*S3  ( P  > 

UV3  <  JQ  >= UV3  <  JO  HU3  <  2  >*S3  <  P  > 

530  UZ3  (  JQ )=UZ3 ( JQ  HUS ( 3  >*S3 < P ) 

N.  END  OF  LOOP  OUER  OFF  BODY  POINTS 
585  P=P+1 
K=K+ 1 
J=J+20 

IF<K-NP  >296,296.509 
P.  END  OF  LOOP  OUER  QUADS 
599  H=AP< I RK  >*DT 

DO  730  I  =  1 , NOBP 
63  FORMAT <2X, I3,3F12.5,9X,4F12.5> 

DX< I >=-(UX I *UX 1 ( I >+UY I *UX2< I HUZ I *UX3< I > > 

DV  ( I  >=-(UX  I  *UV  1  <  I  HUY  I  *UV2<  I  >+U2 1  *UV3< !  >  > 

730  DZ<  l  )=-<UXI*UZU  I  )+>JY  I  *UZ2<  I  HUZl*UZ3<  I  >> 

IF','  IRK.EQ.5)  GO  TO  900 
I F ( I RK . EQ . 4  >  Gu  TO  800 
00  750  1*1. NOBP 
XP<  I  )=XT  ( I  HDX<  I  )*H 
YPC I  )=VT  C I  HDY<  I  )*H 
2P<  I  >=ZT<  I  HD2<  I  >+H 
SKX< I  )=SKXX  !  >+GM< IRK>*DX< I > 

SKY  < I >=$KV< I HOIK IRK  )*0V< I ) 

750  SKZC !  >=SKZ<  I  HGMORK  >*DZ<  I  > 

IRK  =  IRK  +  1 
GO  TO  98 
SCO  H=DT 

DO  S30  1=1, NOBP 

QX<  I  >*-<UX  I  *UX  1  <  f  HUY  I  *UX2<  I  HUZ I  *UX3<  I  >  > 

0V<  !  >=-<UX  i  *UV  1  <  I  HUY  I  *VV2<  I  HUZ I  *UY3<  I  >  > 

0Z< I >=-<UXI*UZU I HUVI*UZ2< I HUZI*UZ3< I > > 

XF‘(  |  >=XT(  I  HH*DX(  I  >/6 .  +SKXf,  I  )*H 
XT ( I )=XF( I > 

VP<  I  >=VT <  I  HH*'DV<  I  >/€• .  ♦SKY(  I  >*H 
YT C : >=VPCI > 

2P( I  >=ZT  ( I  HH*DZ<  I  >/6 .  +SK2C I  HH 
2T< I >=2PC I > 

S'KXX I  >=0. 

SKY ( I  >  =  0. 

S30  SKZ ( I >  =  0 
IRK  =  5 
GO  TO  98 
900  IRK  =  1 

DO  905  1=1, NOBP 

D30=DX( I  )**2+DV< I >**2+DZ< I >**2 
C=\ 1  >=1  -DSO/USQ 
905  CONTINUE 

WRITE '6, 61 )  ITC 

6 1  FORMAT < 6H0  STEF , 14/ > 

WR I TE  <  6 , 62  > 

62  FORMAT C3X.4HL I NE.5X, 1HX, T IX, *HV, T IX, 1HZ,20X,2HUX 

1  10X.2HUZ,  1QX,2HCF  ) 

WRITE <6, 63 >  <1  .XPC  I  >,YPC  I  >,2P<I  >,DX('I  >,DYCI  >.D2' 
WRITE ( 16)  (  XP < I  ) , VP ( !  > , ZP C I > ,  1  =  1, NOBP  ) 

IF< ITC.EQ.NST)  GO  TO  910 
I TC= I TC+ 1 
GO  TO  599 
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910  i F ( i END . EQ . 0  )  GO  TO  23 
REWIND  03 
REWIND  04 
ENDFILE  16 
REWIND  16 
STOP  6 
END 

FUNCTION  SQ2F (X 1 , X2 , Y 1 , V2 Z 1 , Z2 > 
X=X  1-X2 
V=V 1-V2 
2=Z  1-22 

RS=Z**2+Y**2+X**2 

R=RBS <  X X+RBS < V  X+ftBS <  2  >+  1. 0E-20 

R= . 3422*  < R+ ( RS+RS  >/R  > 

R=  R+RS/R 

c,ri7p=  2P^:P.+RP  /p 

RETURN 


APPENDIX  Ull  -  XYZPF  SECTION  PF7 


PROGRAM  PFP7(TAPE7, I NPUT= 128, 0UTPUT= 128, TAPE5= I NPUT, TAPE  17, 

1 TAPE6=0UTPUT , TRPE03 , TAPE04 , TRPE3=TAPE03 , TRPE4=T APE04 ) 

XYZ  POTENTIAL  FLOW  PROGRAM  UERSIGN  4  BHD  UERSION  5  SECTION  7 
COMPUTES  UELOCITIES  AND  PRESSURE  COEFFICIENTS  FOR 
ON  BODY  STRERMLINES 

COMMON  X ( 658 ) , Y  <  658  > , Z ( 658  > , XN  < 650  > , VN < 650  > 

1 ,  ZN  < 650 ) , UX 1 <  650  ) , UV  1  (  650 ) , UZ 1(650), UX2 ( 650 > , UY2( 650 ) 

2,  UZ2C650 ), UX3 < 650 ) , UV3  <  650  > , UZ3C650 ), XC f (658  ), VC  1 (658  > 

3 ,  XC2 <  658 ) , YC2 ( 658 > , XC3 ( 658 ) , XC4 < 658 ) . YC4 < 658 ) , 

4X3  < 653 ) , V3 < 658 ) . Z3<  658  ) , X4 ( 658 ) , Y4  < 658 ).Z4< 658  > 

D I  MENS  ION  XL  ( 1 50  > ,  YL  <  1 50  > ,  ZL  <  1 50  ) .  IJX  ( 1 50  > ,  UV  (  1 50  > ,  UZ  <  1 50  ) , 

1  CP ( 130) , GK  K 1 30  ) , GK2 ( 150  > , H2 ( 130  ) , STHL ( 150  ) , URSS  <150 ) , NOURO <!50> 
5 ,  DMX  ( 650  > ,  F'RuB  <  15  > ,  YC3  ( 658 ) ,  SF  ( 5 ) ,  XCF;  ( 5 ) ,  VCR  ( 5 )  , 

7NSP < 50 > , US  <220  >  , XST ( 50  > , YST ( 50 > , ZST  ( 50  > 

ECU  I URLENCE  (US (20 1 ), HP  >, (VC3, YC2 ) 

RERD(03)(PR0B< i ). 1=1, 15) 

READ ( 03  > ( US ( I ) , I  =  1 , 220  > 

READ (03  )<X< I >,V<  I ).Z< I >.XN( I >,VN< i >.ZN( I ), 

1SKIP, 1=1, NP) 

RER0(C3)SKIP 
lUEFi  =US(220  > 

WRITE (6, 5)  IUER 

5  FORMAT (46H0XVZ  POTENTIAL  FLOW  PROGRAM  SECTION  7,  UERSION  , ;2) 

IF  (IUER.EQ.5)  READ (03)  SLIP 

READ ( 03 )$K I P 

F;EAD(03)SKIF' 

READ(03)SKIP 

RE  AD  ( 03 )  ( UX  1(1),  UV  K  I ),  UZ  1  ( I  > ,  I  =  1 ,  NP  > 

READ ( 03 ) ( UX2 ( ! ) , UY2 ( I > . UZ2 ( i  ) . ! =  1 , NP  ) 

READ (03  >(UX.3(  I  >,  W3(  I  >,  UZ3C  f  >,  I  =  1 ,  NP  > 

REW I  HD  03 
NB=(NP+ 1 1 )/ 12 
DO  80  1  =  1,  HE: 

|RN=I*12 
I S= 1 FN- 1 1 

READ ( C4 )  0,  (XT 1  (..'), VC  1  ( J), XC2( J ), 

1 YC2 ( J ) , XC3 ( J ) , XC4 ( J ) , VC4 ( J ) , X3 (J), V3 (J), 

2Z3  <  J ) ,  X4  <  J  > ,  V4  <  .J),Z4<  J  > ,  <  SK  IP,K=1 ,?>,  J«  IS,  I  FN  > 

NQ=0 

I F ( NQ  .  Nt . IS)  GO  T 0  450 
80  CONT I  HUE 
REw ! ND  04 
DO  on  i=i  nr 

D 1  =  ( XC  H I  )**2+YC 1 ( I  )**2  )* 1.01 
02=  ( XC'2 ( I  )**2  +VC2C I  )**2)*  1.01 
03=(XC3<  I  )*'*2+VC3(  I  y**2  >*  1  01 
04=  ( XC  4 i  )*  ♦  2 ♦VC  4  ( i  )♦  ♦ 2  )*  1  0 1 
90  OMX ( I  )=RMAX  H0 1 , 02 , 03, D4 ) 

1 1  FORMAT  < 3F 12.4.314.F12.4) 

12  FORMAT 12.4,  14) 

M 1 0=75 

1 00  READ (5.11)  UX I , UV I , UZ t.NLIN, MAX J . I  UR  I TE . AMACH 
IF  (EOF'. 5). HE  0  >  NUN=0 


MXJ*MAXJ 

IF  (  MAX..'  LE  0  .OR.  MAXJ 
M I  NJ=M  I  O-MAX'J 
MRX J=M I D+MRXJ 

IF  (MAXJ  GT.MIB+2)  MAXJ=MI0+2 


MAXJ  GT  NF/2)  MAXJ  =  NP/2 


V  V  — 


IF  (  MINU  .LT. 1)  MINJ  =  1 

UR  I TE(7  >  NLIU 

UR  I  TEC  17)  NL I N, UX I , UV I ,  UZ  I 

IF(NLIN.LE.O)  GO  TO  550 

UF;  ITEC6 , 30  )  (PROBC I  ),  1  =  1,  15) 

UFi  I TE ( 6 . 34 )  UXI , UV I , UZ I , NL  IN , MX J ,  I  UR  I TE , AMACH 
34  FORMAT C34H00N  BODV  STREAMLINES  -  INPUT  DATA  /6H  UXI  =,F10.5/ 
16H  UVI  =,F 10.5/6H  UZI  =.F10.5/6H  NLIN=, 110/ 

2  6H  JMAX= ,  1  10. /,SH  IUR1TE=,  I  10,  /,9H  MACH  N0=,F10.5) 
URITE<6,38) 

33  FORMAT  C27H0STREAML I  ME  STARTING  PO I NTS/5H  LINE,  1 IX,  1HX, 52X,  1HV, 
1  12X,  1HZ,  10)  ,  3HNSF' ) 

L I N=NL I N 
DO  45  1=1, LIN 

F;EAD<5 , 12 )  XST <  I  )  ,  VST  Cl),  ZST  Cl),  NSF  c  I  ) 

IF  CE0FC5 ) . EO . 0 .  )  GO  TO  45 
NL I N= I  -  1 

UR  I  TEC 6, 42  )  NL ! N , L I N 

42  FORMAT < 1H0 , i 5 , 2SH  STREAML I NEG  SPEC  I F I  ED  NOT  ,13) 

IFCNLIN.LE  0)  GO  TO  550 
GO  TO  48 

45  UFI  TEC  6. 45)  !  .XSTc 1  >,VST-.  I  >,Z?’-  i  > 

45  FORMAT (IX , 1 3, 2X , 3F 13  5 , 1 9 ) 

48  CONTINUE 

USO=UX I **2+UV I **2+','2 1 **  2 
IF  CAMACH  EO  0.)  GO  TO  1130 
C  COMPUTE  CRITICAL  MACH  NO. 

USD  =  0 

do  iioo  !=i. up 

US  =  (UXi+UXK  I  >+uv|*,..,)2< 1  HUZ'^L'XSC* 

1  <  UX.  I  *UV  1  <  I  >+UV  I  *UV2  <  I  )+UZ  I  *W3  <  I  >  )**2  + 

2  OJX I  *VZ  * ' !  ’  *VZ2  ‘ '  'Z '  * 1  ■>  >•* * 

IF  (US  GT.  USD)  USD  =  US 

1100  CONTINUE 

U  =  SORT  <  USD  /USC*  > 

CMNA  =  1./U 
DO  1110  1=1,3 

CMNE  =  ( ( <CMNA**2+5 . )/6. )**1  75>/U 
CMN-O  =  (((CMNE’ *2+5.  )/6  >**1.75>/U 
IMS  CMNA  =  (CMNR4CMNC-CMNE.'*':*2)/(CMHR+CNNC-2  *'CMME') 

UP.iTEC*  ,  1120)  CMNA 

1120  FORMAT (2 1H  CRITICAL  MACH  NO.  =,F5.3) 

1130  CONTINUE 

C  START  LOOP  ONER  STREAMLINES 
DO'  400  LL=\,NL!N 
0 1 RT=  ? 

101  U I = 1 
RF=1 

'j)  (M I D  )=0 
I.IV<MID)=0. 

UZ  ( M 1 0  )=G . 

CP ( M I D  )=0 
H2'.M !  D  ?=  1 

&  1  ( M ,  D  )=C  . 

GK2 ( M I D )=0 . 

STML''M!D)=0 

102  NQ=NSF(LL) 

LNO=NO 

XL c.M  I  D  )=XST (LL  ) 

VL CM i  0  )=VST ('LL  ) 

ZL (MID  )=ZST < LL  ) 


SEPARATE  CALCULATION  OF  SECOND 
POINT  FRON  MAIN  LOOP 

XLT=(XL(J>-X(NQ)HX3(NQ  H(VLCJ>-VCNQ>>*V3CNQ> 

1  +  (  ZL  (  J  HZ  (  NQ  >  >+23  ( NO  > 

YLT=(XL( J  HXCNQ  >  >*X4(NQ  H(YL<  J  HYCNQ )  >*Y4(NQ ) 

1  +  <ZL<  J  HZCNQ  >  >*Z4(N0 ) 

XLC  J  >=XLT*X3(NQ  HYLT*X4(NQ  HXCNQ  > 

VL  •:  J  )=XLT*Y3  C  m H VLT+V4  (  NQ  >+ V  <  NQ ) 

2L< J  >=XLT*Z3(NQ HYLT*Z4CNQ HZCNQ ) 
i  |qt=moo(NO(4>  +  1 
GO  TO  (630.600,610,620)  IOT 
)  NR=NC!+1 
N  .«=••:  >2 
GO  TO  107 
!  NR=NQ+2 
NU=NQ-  1 
GO  TO  107 
)  NR=N0-2 
NU=NO+ 1 
GO  TO  107 
)  NF;=NQ-1 
NU=NC!-2 

?  UXQ=-(UX ! *0X  1  (NO >+UV I *0X2 (NO HOZ I *0X3 (NQ > ) 

UVO=- < UX I *UV 1 (NO )+0V i *0V2(NQ HOZ I *UY3(NQ > ) 

UZQ*-(OX I *02 1 (NQ >+0V I *022 (NQ HOZ I *023 (NQ > > 

UXF.  -  -  ( OX  I  *0X  1  <  NR  >+UV !  *0X2  ( NR  HOZ I  *0X3  (  NF;  >  ) 

UVR=-(OX  I  *UV  1  (NR  >+UV  I  *0V2(NFl  )+'JZ  I  *UV3CNP, ) ) 

U2F;=-(0X  I  *02 1  (NF;  HOV  I  *022 (NR  HOZ  I  *023(NR ) ) 

UXU=-(OX I *0X 1 (NU )+UV I *0X2 (NU HOZ I *0X3 (NU > ) 

UVU=-(OX I *0V 1 (NU HOV I *0V2(NU HOZ I *0V3(NU > ) 

OZU=- r0X ! *0Z 1 (NO HOV I *0Z2(NU HOZ ! *0Z3(NU ) ) 
TRANSFORM  UELOCITIES  TO  QUAD  SVSTEN 
UQ=0XQ*X3(NQHUVQ*V3(NQHU2Q*Z3(NQ> 

00=UX0*X4 (NQ HUVQ+V4 (NQ HUZQ*Z4(NQ ) 

CSR=  1 . /(XNCNQ  )*XN(NR HVN(NQ HVNCNF; HZN(NQ )*ZN(NR ) ) 
UT=UXF;*X3(NF;  HUVR*V3(NR  HUZR*Z3(NR  ) 

0T= ( UXP*X4 < NR >+UVR*V4 ( NR HUZR+24 ( NR  > >*CSR 
XXR=  <X3(NR >*X3(NQ HV3(NR )*V3(NQ HZ3(NR >*Z3(N0 ) ) 

XVFi=  ( X4  <  NR  )*X3  ( NO  HV4  ( NR  >*V3  ( NO  HZ4  ( NR  )*Z3  ( NO ) ) 

UF; = U  T  *  X  X  R +07  *  X  VR 

VXR*  <X3<  HR  >*X4  ( NQ  HV3  ( NR  >*V4  <  NO  HZ3  <  NP  >*Z4  ( NO  >  > 
V  VF,  -  ( X4  ( NR  >*X4<  HQ  HV4  ('  NR  )*Y4<  HQ  HZ4  ( HR  )*Z4  ( NQ  )  > 
OR=UT*YXR+"T*WR 

0u=UXU*X3 ( NO HUVU*V3(N0 >+020*23 (NQ ) 

CSU=  (XNCNQ  >*XN(NO  HYNCNQ  >*VN(NU  HZHCNG  >*ZN(NU  >  > 

00* ( UX0*X4(NQ >+UVU*V4 ( NO HUZU*24(NQ > )/CSU 
FIND  RELRTIUE  COORDINATES  OF  NEIGHBORING  CURDS 
XD*X  ( NF;  >-X<  NQ  > 

\'D='  '•  NR  >-V(NO  > 

ZD=Z(,;R  >-2(NQ  > 

XT  =X0*X3(NR HVD*V3(NR )+ZD*Z3(NR > 

VTT=XD*X4  (NF;  ,j+YD*Y4(NPi  H2D*Z4(NR  > 

ZT  =XD*XN(NR  HVD*YN(NR  >+ZD*2N(NF;  > 

VT= ( -4*S0RT ( VTT**2+ZT**2 >+VTT*CSR+VTT >*CSR+ . 1666666 
X  p = x T*XXR+YT*X V R 
VF;=XT +VXR+VT  *VYR 
XO=X(HU  HXCNQ ) 

VD=Y(H0>-V(H3> 

ZD=Z(NO  )-Z(NQ > 

XU='(D+X3(N0  >+YD*V3(NQ  H2D*Z3(NQ  > 

VT =>'  Z|+X4(NQ  )+YD*V4(  NO  HZD*Z4(NQ  > 


ZT=XD*XN<NQ )+YO*YN(NU HZD+ZNCNu ) 

VU= ( 4 . *SQRT  < YT**2+ZT**2 )+VT /CSU+YT  >* . 16666667 
FIND  COEFFICIENTS  OF  VELOCITY  FUNCTIONS 
DEN= 1 . /<XR*YU-XU*YR ) 

U 1  =  ( <  UFi-UO  >*YU-  ( U'J-UO  )*VFl )  *DEN 
U2=-( (UFHJQ  >*XU-<UU-UQ )*XR  >*DEN 
V 1 = <  < UR-UQ  >*VU- <  UU-UO )*YR >  *DEN 
U2=-  <  <  UR- VQ  >*XU-  ( UU-UQ  >*XR  >*DEN 
FIND  UELQC I  TV  RT  STREAMLINE  POINT 
USL=UQ+U 1*XLT+U2*YLT 
VSL=UQ+U1*XLT+V2*YLT 
UXP=USL*X3<N0  )+USL*X4(NC ) 

UVP=USL*Y3 < NO >+VSL*V4 ( NQ  > 

L'ZP=!JSL*Z3(NCl  )+VSL*Z4(NQ  > 

FIND  GEODESIC  CURVATURES  GK 1  GK2 
VS0D=USL**2+USL**2 
DEN=USQD+SORT(USQD ) 

GK  1  P=  ( USL4  ‘ VSL*U2-USL*V2  WSl*<VSL*'J  1-USl*V  1  >  )/OEN 
GK2F=CU$L*CUSL"U  1+VSL+U2  (USL^V  J+VSL*V2 )  3/DEN 

FIND  LOCAL  STREAM  FUNCTION 

VI  I.  .111*1  I.".'*  *.  »  -.*■!  \  I - 

V7=U2-UQ:+\'04  r;  1+V2  )/US0D 
XX=U2-CVV-U 1 

C=XLT4U0-YLT*UC-CXV*XLT*YLT-CVV+yLT"’»2-CXX*XLT**2 
:ri2  STREAM  FUNCTION  AT  CORNER  POINTS 


XCR ( 1  >=XC 1 (NO > 

XCR  (2  )=XC2'iN0  .) 

XCP.<3)=XC3<NQ) 

XCR 14  )aXC4<N0 > 

XCR (5  )=XCR( 1 > 

VCR  < 1 >=VC 1 < NO  > 

VCR  <2  '« 

VCR  (3  l)=VC3(N0 ) 

VCR  (4  )=VC4lfiQ ) 

VCF;<5  )=VCR<  1 ) 

DO  110  N=t,4 

'  SF  IN  >=CO-VO*XCFliN  >+UO*YCR<N  >+CX V*XCP < N  >*  VCR < N  >+CVV*VCF;(N  >** 
1  +CXX*XCR<N>**2 
Srs5>SF<  1  > 


DC  1 

20  N=  1 , 4 

if  <; 

SF  <N >*SF CN+  1  > 

GE.  0.  >  ( 

30  TO  120 

YM-  ,• 

XcF,  ;N  .)+xlR'-.N*  1  > 

V*:  =1 

V(1=c 

VCR  i  N  >+l.CR(N+  1  i 

)* .  5 

F  i  NO 

INTERSECTION  W 

ITH  SIDE  Oc 

QUAD . 

L'pMr 

C  0  -  V0-f  X  •  ■  !j  C  *  \  7  H  C  XV*  XN  *  V?'1 ♦ C ; 

r|y*VM**2+; 

AC =2 . * ( SF < N  >-2 . "SFN+SF <N+ 1 > > 

E:C= .  *$F<N  >-4 .  *SFM+SF  <N+ 1 ) 

IF  IRC  .EC.  0)  GO  TO  112 
^  *fcf>.=  r ,  n  >  y 

TP=(B'C+SR >/<2 . "AC ) 

IF  <TF  LE  1  AND  TP  GE  0  >  GO  T0  115 

TP=f  E'C-S0  '  /fl2  *AC  ) 

GO  TO  115 

113  IF  <BC  EO.  0  >  GO  TO  120 

TP=SF<N)/E:C 

1 15  XHP=(  1 .  -TP  )*XCRtN  )+TP*XCF;CN+ 1 ) 

VHP® (.  1 .  -TP  >*YCR<N  >+TP*YCR(N+  1  ) 

TESTP=( (XHP-XLT )*UO+<YNP-VLT >+00 )*D I RT 
IF  (TESTF  LE  TEST)  GO  TO  120 
test=testp 
XNT =XHP 


YNT=YNF 
120  CONTINUE 

IF<  TEST  EO.  CO  GO  TO  280 

fi'JERAGE  LAST  UELOCITV  AND  CUR'JATURE 

UX<  J  >=CUX< J  >+UXP  >*AC 

UV ( J  )=  C UV ( J J+UYP >+hF 

UZ  C  J  >*  C  UZ  C  J >+UZP >*AF 

OK K J  )= <  OK  1 ( J )+GK 1 F  >*FIF 

GK2C J  >=CGK2C J  >+GK2P )*AF 

H2< J >=H2< JL >*<2. -OK  1 <JL >*<STHL<J )-STML<UL  )> >/<2 . +GK 1  ( J >+ 
1  STML  C J )-STNL C JL )  > > 

CPC J  )=  1 . -CUXCJ  )++2+UVC J  )**2+UZC J )++2 )/USO 
DABS  C  J  )=SQFlT  C 1 .  -CPC  J  >  > 

COMPUTE  UELOCITV  AT  NEXT  POINT 
NQUAD C J >=N0 

JL=J 

J=J+JI 

USL=UQ+XNT+U 1  +YNT*l»2 
USL=UQ+XNT*U !+VNT*V2 
UXCJ  )=USL*X3  C  NO  )+USi->%4CNQ  > 

I  j1  lc,  j  ';=(  !«•!  *U'.;  f  Mf,  i+i  IC;I  *y4  (W%  s 

UZ  C  U  )=USi_*Z3  C  NO.  >+U  S_:fZ4  C  NO  > 
uuHPUTE  GEuufc  i  CUAUhTUhEc. 

UFOr!=I  !'-'L**2+USL*!*:2 
OEN=USQD*SQF;T  CUSOD  ) 

OK  1<J  >=  usl* C USL*U2-USL*U2  >-U$L*CUSL*U  1— US  l*U  «  >  >/PEw 
GK2C  J  }= ( USL* < USL*”  *  ♦VS  L*  U2  ■"  L'SL  *0 1  ■‘•US  L  *v2  >  )/CEN 

COFiO=SORT C  C  XNT-X.LT  >*+ 2+  C VNT-VLT  )**2  >*D  i  FiT 
STML  <  J  >=3TML  < JL  >+C0F;0 
C0“1=-.’E  H2 

H2 <  J  >=H2 C  JL  >* < 2 .  -CORD*GK  1 C  JL  >  > /  C 2  .  +COF;D*GK  1 C  J  > ) 

(*r,  -  »  <C  /)  ic  “> 

LiREliC  J  )=SQF:T  C 1  -CF  C  J  )  > 

AF=  5 
LNO=NQ 

XL  C  J  >=X  NT'fXS  i  NC|  >+VNT*  X'4  C  NO  )+X  C  NO  > 

Vj.it  j  )sXMT*S'3CN0  >+YNT*Y4CN0  >+VCNO  ) 

ZL  C  J  )=XNT*Z3  C  NO  )+VNT+Z4  C  NO  )+Z  <  NO  > 

!P  C  J  .  LE .  NINJ  OF,  J  GE  UAXJ  )  GO  Tu  280 
FIND  NEXT  QUAD. 

1=1 
i50  NO* I 

i PC!  EO  LNC >  GO  TC  280 

TEST=CXLCJ  )-?(  I  ) >**2+<Yu<J  >-V<  l  )  >**2+ 

1<ZLCJ>-2c  I  > )*+2-ZW '  1  ) 
iFCTES”  jT . 0 .  )  GO  TO  230 
DS  1=\XC  t  C  i  )-XC2C  I  >  >*+2+CVC  I C  i  )-VC2\  i  )  >*+2 
DE2=CXC2C  I  )-XC3f.  I  )>**2+(VC2'‘  !  )-723c :  > 2 
OSS- CXCS C  i  )-XC4C  I  ))**2+<VC3(  I  >-VC4C  I  )  )*  ♦  2 
DS4*<XC4<  I  >— XC 1  < !  >>**2+<VC4<  l  >— S-'C  1  ■.  i  >  >+*2 
':>  4_T= '"’XL  ( J  V-Xc  1  )  )*\3(  J  >+  C  VL  ( J  >-V <  J  )  i 
1CZLC J >-ZC I  )  )*Z3C I  > 

\>lt = ( X'L r-  J  ?—  x  (  |  i  r  4 :  t  ’?+  ( VL  (  J  )-V f,  I 
1CZL(J>-ZC i )X*Z4;  i > 

ZLT*<XL<J)-X< I  >>*XN< I HCYLC J )-Y< I >)*VfK  I H 

KZ’_cj)-z:  ’  ' 1  ■) 

ZS0=ZLT**2 
TEST=ZSO- . i+DMi' I > 

IFCTEST  GT.O  >  GO  TO  280 

AC  1=SQRT CZSO+CXLT-XC 1  <  I  )  v+’+2+CVLT-VClf.  I  )) 

1**2) 

RC2=SQF:TCZS0+CXLT-XC2C  I  >)*+2+,:YLT-YC2<  >  >>*+2 . 


RC3=SQRT  CZSQ+ (XLT-XC3C I )>**2+<YLT-VC3< ! >)**2> 
RC4=SQRT<ZS0+<XLT-XC4< I >  )*+2+<YLT-YC4< I >>**2> 

TEST=  ( <RC 1+RC2 )**2  >-DS 1  *1.21 
IF(TEST . LT .0 .  )  GO  TO  105 
TEST*  < (RC2+RC3 >**2  >-DS2  *1.21 
IF<TEST.LT.O. )  GO  TO  105 
TEST=  < (RC3+RC4 >**2 >-DS3  *1.21 
IFCTEST.LT.O. )  GO  TO  105 
TEST*  ( (RC4+RC 1 )**2  >-DS4  *1.21 
1F<TEST.LT .0. >  GO  TO  105 
i  1=1+1 

IF< I  LE.NP)  GO  TO  250 
:  IF  (DIRT  . LT.  0.  )  GO  TO  285 
D I RT=- 1 . 

J !  =- 1 
JMRX=J 
GO  TO  102 
I  JM I N= J 
SS=STML(JMIN) 

DO  290  J*  JM  I N ,  JMfi.X 
i  STML ( J  >=STML ( J )-SS 
JMN=JMfN+1 
JMX-JMfiX-2 
RF  =  1 . 

L=JMN 

UR!TE<6,30><PR0B<  I  >,  1=1, 15) 

!  FORMAT <  1H 1 ,  15ft 4 ) 

WRITE  <  6 , 20  >UX l , UV I ,  U2 1 

»  FORMAT <  18H0  ONSET  FLOW,  UX I = ,F0 . 3, 2X. 4HUVI =, F6 . 3, 2X, 4HUZ I = , F6 
UR  ! TE  <  6 , 50 )  LL , MSP (LL>,XST(LI_>. VST ( LL  > , 2ST < LL  > 

)  FORMAT  < 1 1 HO  LI  HE  NO .  ,12,31m  PASS  I NG  THROUGH  QUflDR 1  LATERAL  , I 

1  28H  WITH  STARTING  POINT,  X=,  F 12 . 5 , 2X , 2HV=, F 12 . 5, 2X , 

2  2H2* ,F  12 . 5  //) 

IF  <JMIN.LE.MINJ  .OR.  JMft.T . GE . MAXJ >  WRITE <5  55'* 

•  FORMAT OSH  PROBABLE  ERROR  -  LINE  IS  VERY  LONG) 

DO  330  J=Jf1N,  JMX 

IF  <  <STML< J+2 >-STML(L- 1 > ) . LT .  8.*<STML(J+1 )-STML  <L)>)  GO  TO 
UR  I TE  <6.31 0 )  XL < L ) . VL < L  > , 2L  <  L  > , XL  < J+ 1 ) , VL < J+ 1 ) , ZL <  J+ 1 ) 

•  FORMAT < 14H  PO I  NT  DELETED  ,  10X , 3F  12.5, 10X , 3F 12 . 5  > 

STMLCL )=<RF+STML<L )+STML<J+ 1 ) )/<RF+ 1 . ) 

XL<L  )=<AF*XL<L  )+XL< J+ 1 ) >/<RF+ 1 . ) 

VL < L  >=<FlF*YL(L  >+VL<  J+ 1 )  >/<RF+ 1 .  ) 

ZL<L )=<RP+ZL(L )+2L< J+ 1 > >/<AF+ 1 . > 

UX < L >= ( RF+UX < L >+UX< J+ 1 > )/<AF+ 1 .  ) 

UV <L )=< RF  *UV < L )+UV <  J+ 1 ) > / <RF+ 1 . > 

UZ(L  )=<flF*lJZ<L  >+UZ<  J+ 1 )  >/<A-+ 1 .  > 

OK  KL  )=<AF*GK  1  <L  >+GK  1  <  J+  1  >  >/<AF+ 1 .  > 

GK2CL )= ( RF+3K2 C L )+0K2 ( J+ 1  ))/CAF+ 1 .  ) 

H2 <L)=< AF*H2 < L  >+H2 < J+ 1 ) )/< AF+ 1 .  ) 

CPCL  )=  1 . -<IJX<L  >**2+UV<L  )**2+UZ<L  >**2 )  /USQ 
UABS ( L  >=SQRT < 1 . -CP <L ) > 

AF=RF+ 1 . 

GO  TO  320 
)  RF=  1 . 

L=L+  1 
K=J+ 1 

STML<L)*STML<K> 

XL<L  >=XL<K ) 

VL(L  >=VL<K ) 

ZL<L  )=ZL<K ) 

UX<L  )=UX<K ) 

UV(L  >=UY(K ) 


o  o  o 


UZCL  )=UZ<K ) 

GKKLMSKKfO 

GK2CL)=GK2<iO 

H2CL)=H2(K) 

CP(L)=CPCK) 

UftBSCL)=UA6S(K) 

NQURD  C  L )=NQUflD  <  K  > 

330  CONTINUE 
L=L+ 1 

STML  (  L  >»STML  <  JMflX  > 

XL(L)=XL<JMRX) 

VL(L)=VL(JMRX) 

ZLCL  )=ZL< JMflX  > 

UXCL>=UX<  JMflX) 
uv<L>=uv(jm:o 

UZCL)=UZCjriRX) 

GK  KL  )=GK  1  (UMAX ) 

GK2  <  L  >=GK2  ( JMflX  ) 

H2CL)=H2< JMflX) 

CPCL >=CPC JMA.X' ) 

URBS  C L )=URBS < JMRX  > 

jiiflx=L 

NQURD <  JMflX  >=NQIJflO<  JMflX- 1  ) 

NQURD  ( JM I N  )=NQLIflD  <  JM !  N+ 1 ) 

UR!TEC6,51) 

5 1  FORMAT  < 4H0  1 . 6X . 1HX . 9X , 1HV , 

1 9X  ,1  HZ ,  09X ,  2HUX ,  8X ,  2HU  V ,  8X ,  2HUZ ,  09X , 

22HCP ,  8X , 2HK 1 .  8X, 2HK2 .  8X , 2HK2 , 8X , 2HSL ,8X, 1HU, 9X , 1 HP  > 
!F  CflMftCH  .EO.  0.  )  GOTO  1 160 
C  ***  COMPUTE  COMPRESS ! B I L i TV  CuRREuTluN 
DO  1150  J=JM IN, JMflX 

USD  =  < UX ( J )**2+UV <  J >**2+U2 < J  >**2 > /USQ 

USDft  =  USD 

SM  =  fiMflCH**2 

DO  1140  1=1,3 

P,  =  <  1 .  +  .  2*SM*  < 1 .  -USDR ) ) 

IF  <F;  .LT.  .000001)  R  =  .000001 
USDS  =  USD/R**2.5 
R  =  <  1.+.2*SM*< 1.-USDB)) 

IF  CR  .LT.  .000001)  R  =  .000001 
USOC  =  USD/R**2 . 5 

1140  USDR  =  <USDC*U$Dfl-U$DB**2 )/vUSDC+USDA-2 . +USDB ) 

R  =  < 1 . + ,2*SM*< 1.-USDA)) 

IF  <R  .LT.  .000001)  R  =  .000001 

R  =  1 . 25 

IJXC  J >  =  UX(J )/R 

UVCJ)  =  UV<J)/R 

UZCJ)  =  UZC  J )/R 

UftBSCJ)  =  SORT (USDR) 

1 150  CPC J >  =  <R+*2 .8-1. )/<. 7*SM ) 

1150  CONTINUE 
K.=0 

DO  53  I =JM IN, JMflX 
K=K+  1 

53  UR i TEC 6, 60 )  K.XLCI ) , VL < 1 ),ZLCI ),UX<I ).UVCI >,UZCl >,CPCI  > 
1  GK  1  (  :  ) .  GK2  <  I  ) ,  H2  C  I  ) ,  STML  C I  ) ,  URBS  CD.  NQURD 

60  FORMAT (IX. 1 3, 3F 10 . 5, IX, 3F 10 . 5 , IX, 6F 10.5, 16) 

8  FORMAT <3F 12.5) 

UR  I  TEC  17)  K,  CXLC I ),VL( I ),ZLC I ),NOUHD< I ),  l=JMIN, JMflX 
(WRITE  ,LE .  0  —  WRITE  SL,U,H2,K2 

I  WRITE  GE.  2  —  WRITE  X,V.2,CP 

IURITE  .EQ.  1  --  WRITE  SL,U,H2,K2  AND  X,V,Z,CP 


v  *  -•  *  ~.’fTTTyT>l  * '  -  v  •*.'■."*  i'*V*  l^1 7* r*r iwvi  t*  '.■V'.ws  *v* 


IF  CiUfiiTE.GT. 1 )  Gu  TO  340 

URITEC7)  K,  CSTMLC I >,URBS< I  ),H2< I ),GK2< I ), l=JMIN, JMAX) 
340  IF  < IWR1TE.LT. 1)  GO  TO  400 

UR I  TEC?)  K,  <XL< I  )/VLC I  ),ZLC I  ),CP< I  ),  I = JM I N . JMflX ) 

GO  TO  400 

300  UR  I  TEC 6, 50)  NSPCLL) 

UR  I TEC6, 65  > 

GO  TO  282 
400  CONTINUE 
GO  TO  100 

C  READ  NEXT  SET  OF  STREAMLINES 

450  URITE<6,451 )IS,NQ 

451  FORMftT C 14H  TAPE  04  ERROR, 21 4) 

550  ENDFILE  7 

FOUND  7 
ENDFILE  17 
REUIND  17 
REWIND  04 
STOP  7 
END 


I  30 


flPPEHO 1 X  Ulll  -  TRIRXIRL  ELLIPSOID  INPUT  FILE 

SAMPLE  PROBLEM  TR! AXIAL  ELLIPSOID 

2S0  5  150  150  150  3  .00001  00000 


5  150  150 

150  3 

.00001  0 

0 

0 

0 

0 

0  0  0 

.00000 

. 00000 

.00000 

1 

1 

1 

0 

. 00000 

97861 

. 00000 

.  10286 

1 

2 

1 

0 

.00000 

.  90789 

. 00000 

. 20960 

1 

O 

1 

0 

. 00000 

.  82360 

.00000 

.28359 

1 

4 

1 

0 

. 00000 

.71583 

.00000 

.34914 

1 

5 

1 

0 

. 00000 

. 62478 

. 00000 

. 39040 

1 

6 

1 

0 

. 00000 

. 52537 

. 00000 

. 42544 

4 

1 

7 

1 

0 

. 00000 

.44721 

.00000 

.44721 

1 

8 

1 

0 

. 00000 

.  37530 

. 00000 

. 46345 

1 

9 

1 

0 

. 00000 

.  29822 

. 00000 

.47725 

1 

10 

1 

0 

. 00000 

.23692 

. 00000 

. 48576 

1 

1 1 

1 

0 

. 00000 

. 16957 

. 00000 

.49275 

1 

12 

1 

c 

CHHH3 

. 1 1467 

.00000 

.49670 

1 

13 

1 

0 

.00000 

. 05243 

. 00000 

. 4993 1 

1 

14 

1 

0 

. 00000 

.  00000 

. 00000 

. 50000 

1 

15 

1 

0 

. 00000 

.  99875 

. 10000 

. 00000 

2 

1 

1 

0 

. 00000 

97739 

.  10000 

. 10273 

2 

2 

1 

0 

. 00000 

.  90676 

.  10000 

. 20934 

2 

3 

1 

0 

. 00000 

.  82257 

.  10000 

28323 

2 

4 

1 

0 

. ooooo 

.71494 

.  10000 

. 34870 

ni 

5 

1 

o 

. ooooo 

.62299 

.  IQuOu 

. 3899  1 

o 

6 

1 

0 

.00000 

.  5247 1 

. 10000 

. 42490 

o 

4m 

7 

1 

o 

. ooooo 

44665 

.  10000 

.44655 

2 

8 

1 

0 

. ooooo 

.  o  »•  4  c  o 

.  IG'JU'J 

.46287 

4m 

Q 

4 

l 

o 

.00000 

.  29785 

.  10000 

. 47665 

2 

10 

1 

0 

. ooooo 

.  23663 

i  n~:r:ri 

. 485 16 

2’ 

11 

1 

o 

. 16946 

.  10000 

.49213 

2 

12 

1 

0 

. ooooo 

. ! *453 

.  10000 

.49608 

**. 

4m 

13 

1 

0 

. ooooo 

. 0524  1 

.  10000 

. 49859 

2 

14 

1 

0 

.  ooooo 

.  00000 

.  10000 

.49937 

2 

15 

1 

o 

. ooooo 

. 99459 

, 00000 

3 

1 

1 

0 

. ooooo 

.  9737 1 

. 20000 

.  10234 

3 

2 

1 

0 

. ooooo 

. 90334 

20000 

20855 

3 

3 

1 

0 

.00000 

. S 1547 

. 20000 

.28217 

3 

4 

1 

0 

. ooooo 

.71225 

.  20000 

. 34739 

3 

Ej 

1 

0 

. ooooo 

<e  ^  4 

. 20000 

. 38845 

3 

6 

1 

0 

. ooooo 

.52274 

. 20000 

. 42330 

3 

7 

1 

0 

. OuOuG 

44497 

20000 

.44497 

3 

8 

1 

0 

. ooooo 

.  i-  ( 4  2 

.  2 . vJU 

.  4f.  1 13 

3 

9 

1 

0 

. ooooo 

. 29672 

. 20000 

. 47486 

o 

10 

1 

o 

.  OlAW 

. 22574 

.  20000 

. 48333 

3 

1 1 

1 

1J 

flfl! 

159S2 

. 20000 

.  4yU2!b' 

•“> 

12 

1 

o 

.  UJ.V'J 

.  1 1410 

. 20000 

. 4942 1 

3 

13 

1 

0 

. ooooo 

. 05222 

. 20000 

4958 1 

yjl 

14 

1 

0 

. ooooo 

.  00000 

. 20000 

.49749 

O 

15 

1 

0 

. ooooo 

Q8SF.Q 

. 30000 

, 00000 

4 

1 

1 

0 

. ooooo 

. 96754 

90000 

10169 

4 

i. 

1 

0 

CfSOOO 

. 89762 

. 30000 

•*>  “i  —i  ■•i  •; 
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